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We compute the conservative piece of the gravitational self-force (GSF) acting on a particle of 
mass mi as it moves along an (unstable) circular geodesic orbit between the innermost stable orbit 
and the light ring of a Schwarzschild black hole of mass m,2 3> mi. More precisely, we construct 
the function h^ L (x) = hfy ,u^u v (related to Detweiler's gauge-invariant "redshift" variable), where 
h^i, L (cc mi) is the regularized metric perturbation in the Lorenz gauge, is the four- velocity of 
mi in the background Schwarzschild metric of ni2, and x = [Gc~ 3 (mi + m,2)f2] 2//3 is an invariant 
coordinate constructed from the orbital frequency fi. In particular, we explore the behavior of h^u 
just outside the "light ring" at x = ~ (i.e., r = 3Gm2/c 2 ), where the circular orbit becomes null. 
Using the recently discovered link between and the piece a(u), linear in the symmetric mass 

ratio v = mim2/(mi + m.2) 2 , of the main radial potential A(u,v) = 1 — 2u + va(u) + 0(v 2 ) of 
the effective one body (EOB) formalism, we compute from our GSF data the EOB function a(u) 
over the entire domain < u < | (thereby extending previous results limited to u < |). We 

find that a(u) diverges like a(u) w 0.25(1 — 3u)~ 1//2 at the light-ring limit, u — > (|) , explain the 
physical origin of this divergent behavior, and discuss its consequences for the EOB formalism. We 
construct accurate global analytic fits for a(u), valid on the entire domain < u < | (and possibly 
beyond), and give accurate numerical estimates of the values of a(u) and its first three derivatives 
at the innermost stable circular orbit, as well as the 0(u) shift in the frequency of that orbit. In 
previous work we used GSF data on slightly eccentric orbits to compute a certain linear combination 
of a(u) and its first two derivatives, involving also the 0(v) piece of a second EOB radial potential 
D{u) =1 + 1/ d(u) + 0(u 2 ). Combining these results with our present global analytic representation 
of a (it), we numerically compute d(u) on the interval < u < |. 

I. INTRODUCTION 

For much of its long history, the two-body problem in general relativity has been studied primarily within two 
analytical approximation frameworks, one built around the weak-field limit and the other around the test-particle 
(geodesic) limit. The first analytical framework, formalized in post-Newtonian (PN) and post-Minkowskian theories, 
is (a priori) applicable only when the two components of the two-body system are sufficiently far apart. The second 
analytical framework is (a priori) relevant only when one of the masses is much larger than the other, in which 
case the dynamics can be described, at first approximation, as a geodesic motion on a fixed curved background. 
Recently, rapid developments (mixing theoretical and numerical methods) in the field of gravitational self-force (GSF) 
calculations (see [T] for a review) have allowed one to go one step beyond the geodesic approximation, giving access 
to new information on strong-field dynamics in the extreme-mass-ratio regime. In addition, since 2005 it has been 
possible to accurately describe the coalescence of two black holes of comparable masses by using three-dimensional 
numerical simulations based on the fully nonlinear Einstein equations. The progress in interferometric gravitational- 
wave detectors has brought with it the imminent prospect of observing gravitational radiation from inspiralling and 
coalescing astrophysical binaries, and with it the need to compute, in an efficient and accurate way, the form of 
the many possible gravitational-wave signals emitted by generic binary systems (having arbitrary mass ratios and 
spins, and moving on generic orbits). It has become clear over the past few years that the best way to meet the 
latter theoretical challenge will be to combine knowledge from all available approximation methods: PN theory, 
post-Minkowskian theory, GSF calculations, and full numerical simulations. 

Within this program, the effective one body (EOB) formalism [2]-[5] was proposed as a flexible analytical framework 
for describing the motion and radiation of coalescing binaries over the entire merger process, from the early inspiral, 
right through the eventual plunge and final ringdown (see Ref. [6] for a review). The central posit of the EOB 
formulation is a mapping between the true dynamics and an effective description involving an effective metric, together 
with an extra "mass-shell deformation" phase-space function Q involving (effective) position and momentum variables. 
If the two objects are nonspinning black holes with masses mi and m.2, then in the extreme-mass-ratio limit [i.e., when 
the symmetric mass ratio v = mirri2 / (mi +m2) 2 tends to zero] the effective metric is expected to reduce smoothly to 
the Schwarzschild metric, while Q must vanish. For a general mass ratio, i.e. for a non-zero value of v in the interval 
< v < 1 the effective metric involves two initially unspecified functions of two variables ("EOB potentials"), 
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denoted A(u;v) and D(u;v). Here u is the dimensionless "inter-body gravitational potential" u = GM/(c 2 rEOB), 
where M = m,\ + 777,2 denotes the total mass, and teob is the (EOB-defined) radial separation between the two 
objects. In the current, "standard" formulation of the EOB formalism, the motion in strictly circular binaries is 
governed by the potential A(u, v) alone. The (conservative) dynamics of slightly eccentric binaries involves, besides 
A(u,v), the second EOB potential D{u;v). More generally, the conservative dynamics of arbitrary orbits (described 
by the full EOB Hamiltonian) involves, besides A(u; v) and D(u] v), the third EOB function Q(u,p v ,p r ; v), which a 
priori depends on the four variables (u,p v ,p r , v), where p v is the angular momentum and p r is the radial momentum 
canonically conjugated to the radial variable teob- 

Post-Newtonian theory only gives access to the expansions of the EOB potentials in powers of the inter-body 
gravitational potential u, while keeping the exact dependence upon v. For instance, PN calculations at the third PN 
(3PN) approximation lead to the exact knowledge of the coefficients A 2 (y),A^{u) and A^(y) in A(u, v) = 1 — 2u + 
A2(v)u 2 + A^(v)u 3 + A^v)^ + 0(u 5 \nu) , with the remarkably simple result [3] that the 1PN coefficient A^iv) 
vanishes, and that the 2PN, A3 (v), and 3PN, Ai(y), coefficients are both linear in v [thanks to some remarkable 
cancellations; the function A^iv), e.g., is a priori a cubic polynomial in v\. In order to apply the EOB formalism 
to the description of the final stages of coalescing binaries, it is necessary to somehow improve the behavior of these 
(weak-field; u <C 1) PN expansions, and to extend the knowledge of the functions A(u), D(u) into the strong-field 
regime u = 0(1). Two different methods have been proposed to perform such a strong-field extension. Both methods 
exploit the flexibility of the EOB framework, which naturally allows for either the introduction of unknown parameters 
(parametrizing higher-order PN terms), or for the introduction of unknown functions (linked to GSF theory). 

The first method used for "upgrading" the PN expansions of A(u; v) and D(u; v) into functions which are (tenta- 
tively) valid in the strong-field regime 77 = 0(1), was to replace them with suitably resummed expressions, namely 
some Pade approximants of either the currently known PN expansions |4 , or of PN expansions incorporating some 
undetermined coefficients parametrizing as-yet-unknown higher-order PN terms [THlOj. As results from strong-field 
numerical relativity (NR) simulations started to emerge, it became possible to "calibrate" some of these unknown 
parameters, by finding the values that "best fit" the NR data fTHTO]. The resulting NR-fitted EOB formalisms have 
been found to provide a useful analytic approach to the two-body problem in both the weak- and strong-field regimes 
and across all mass ratios [9TTT5]. 

The second method for extending the validity of the PN expansions of A(u; v) and D(u; v) is to use information 
from GSF theory [16]. Essentially, while PN theory (in the EOB context) involves the expansion of A(u; v), D(u; v) 
and Q(u, p^jPr', v) in powers of u (for fixed v), GSF theory involves the expansion of these functions in powers of v (for 
fixed it). For instance, the GSF expansion of the A potential is of the form A[u; v) — I ~2 u + v a(u) + v 2 a2(u) + 0(i/ 3 ), 
while that of D(u; v) starts as D(u; v) = 1 + v d(u) + v 2 c?2(w) + 0(v 3 ), where we suppressed, for notational simplicity, 
the index 1 on the coefficients a(u) and d(u) of the first power of v ("first GSF level"). Note that all the GSF 
coefficients a(u) , 02(74) , d(u) , d<2(u) are functions of u, and are a priori defined for arbitrary values of u, including 
strong-field values u = 0(1). Since 2008, calculations of the GSF in Schwarzschild geometry are providing valuable 
information on various invariant aspects of the post-gcodesic dynamics in binaries of extreme mass-ratios. This offers 
a new opportunity for improving the EOB formalism by acquiring knowledge on the strong-field behaviour of the 
various functions a(u), 02(77), d(u), d 2 (Tt), ... . The GSF data are particularly useful for this purpose since they are 
highly accurate (GSF calculations involve only linear differential equations), and because they give access to a portion 
of the parameter space inaccessible to either PN or NR: strong-field inspirals in the extreme-mass-ratio domain. 
Furthermore, in GSF calculations (unlike in NR) it is straightforward to extract the conservative (time-symmetric) 
aspects of the dynamics separately from the dissipative ones. This is an advantage because the two aspects are dealt 
with separately in EOB. 

The promise of such a GSF-improved EOB formalism was first highlighted in Ref. [16] . That work suggested several 
concrete gauge-invariant quantities characterizing the conservative dynamics of the binary, which can be constructed 
(in principle) using knowledge of the GSF, and would provide accurate information about the ^-linear EOB functions 
0(77), d(u). As a first example, Ref. [H] used the GSF computation [T7] of the 0(v) shift in the value of the frequency of 
the innermost stable circular orbit (ISCO) of the Schwarzschild black hole, to determine the value of the combination 
a(u) = a(u) + ua'(u) + ^7t(l — 2t*) a" (77) (where a prime denotes d/du) at the ISCO potential value u = |. Ref. [T5] 
also proposed that a GSF computation of the frequency and angular momentum of a marginally bound zoom-whirl 
orbit could be used to determine the separate values of 0(77) and a' (77) at the much stronger-field point u = \ (the 
"whirl" radius), but such a computation is yet to be performed. 

More importantly, Ref. [11] has shown that a computation of the GSF-induced correction to the periastron advance 
of slightly eccentric orbits along the one parameter sequence of circular orbits, would allow one to compute the 
combination |51j p(u) = a(u) + (1 — 677) g?(t7) as a function of u over the entire range where circular orbits exist, i.e. 
< 77 < |. The calculation of the EOB function p(u) was then performed in Ref. [18] along the sequence of stable 
circular orbits (i.e. < 77 < g), when computational tools for the GSF in eccentric binaries became available [19] , 
Ref. [18] also made the following important point. By combining PN information about the behaviour near u = 
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of functions such as a(u) or p(u), together with the GSF-computed values of these functions at a (possibly sparse) 
sample of strong-field points u = ui, U2, ■ ■ ., one can construct simple (Pade-like) analytic representations which can 
provide accurate global fits for the corresponding EOB functions. Then, in turn, these global representations can be 
used to analytically represent other GSF functions of direct dynamical significance. Ref. [T5] demonstrated this idea 
by constructing a simple, yet accurate, global analytic model for the periastron advance in slightly eccentric orbits, 
using only a small set of strong-field GSF data in conjunction with available weak-field PN information. In subsequent 
work |14j this model was successfully tested against results from fully nonlinear numerical simulations of inspiralling 
binaries. 

Unfortunately, knowledge of the GSF- induced periastron advance only gives access to the combination p{u) involving 
the functions a(u), a'(it), a"{u) and d(u), and it is not sufficient for determining the individual potentials a(it) and 
d(u) separately. This situation was cured in recent work by Le Tiec and collaborators. In Ref. [3D] Le Tiec et al. have 
"derived" (using a mixture of plausible arguments) a "first law of binary black hole mechanics" , relating infinitesimal 
variations of the total energy £ and angular momentum J of the binary system to variations of the individual black 
hole "rest masses", and (for mi <C ttiq ) involving Detweiler's red-shift variable z\ associated with m,\ [21]. (The 
validity of this relation was established rigorously only through 3PN order.) Based on this relation, further work 
[2"2l about the functional link between £, J and the dimensionless orbital frequency parameter x = (GMfl/c 3 ) 2 ^ 3 led 
Barausse, Buonanno and Le Tiec [23] to derive a simple direct relation between the 0{v) piece of the function z\{x) 1 
and the 0(v) EOB function a(u) (evaluated for the argument u = x). This relation shows that GSF calculations of 
the 0{v) piece of the redshift function z\(x) of mi, along circular orbits, allows one to compute the function a(u), 
separately from the second 0{v) EOB function d{u). [Using the (quite simple) EOB theory of circular orbits, it is 
then easy to derive from a{u) the functions relating £ and J to both u and x; see Refs. [TBI [23] an d below.] By putting 
together the so-acquired knowledge of the function a (it) with Ref. [TB]'s GSF computation of the combination p(u), one 
then has separate access to the second EOB function d(u), thereby completing the project initiated in Refs. [TBI, ITS] 
of using GSF data to determine the (separate) strong-field behaviors of the two main 0(y) EOB potentials a(u) and 
d(u). Note, however, that this still leaves out the third EOB function Q(ii, p^p,.; v). 

The analyses of Refs. [22] and [23] relied on numerical GSF data for zi(x), which have so far been available only 
for x < 1/5 [3TJ[53]. This allowed the determination of the EOB potentials (and of £ and J) through 0(v) only in 
the restricted domain < u < 1/5. The EOB potentials remained undetermined in the strong-field domain u > 
In the extreme-mass-ratio case, this domain corresponds to the region r < hGm-ijc 2 outside the large black hole 
of mass 7712 (where r is the Schwarzschild radial coordinate associated with m,2, which coincides with teob in the 
mi — > limit). Note that the gravitational potential varies steeply in this region, so that the EOB functions might 
well vary correspondingly fast and possibly in a non-trivial way, potentially giving rise to interesting new physics. In 
this regard, we emphasize that (as is clear in EOB theory) it is the gravitational-potential coordinate u, and not r 
itself, which best parametrizes the strength of the gravitational field. We note in this respect that the gravitational 
potential difference across the seemingly "small" domain extending between the ISCO and the light ring (below which 
there exist no circular geodesic orbits), 3Gm 2 /c 2 < r < 6Gm 2 /c 2 (i.e., | < u < |), is as large as that across the 
entire domain 6Gm 2 /c 2 < r < 00 (i.e., < u < g). This lends a strong motivation for extending the analyses of Refs. 
[2"2] and [23] to the domain | < u < §. 

In this paper we obtain numerical GSF data for z\(x) for circular geodesic orbits with radii in the range 3Gm 2 /c 2 < 
r < 150Gm 2 /c 2 . We use these data to compute the numerical values of the function a(u) on a dense set of 11 values, 
extending down to 11 = |. We then construct a global analytic fit for the function a(u), valid uniformly on < u < |. 
We pay particular attention to the behavior near u = |, which, in the limit v — >• 0, represents the light ring (LR) of 
the Schwarzschild black hole, where circular geodesic orbits become null. 

It should be commented immediately that the interpretation of the GSF near the LR is a subtle one: for any 
finite (nonzero) value of v, there are sufficiently small values of u — | for which the mass-energy of the small particle 
becomes comparable to that of the large black hole, at which point perturbation theory breaks down and the GSF 
approximation ceases to be meaningful. In principle, however, it is possible to make the GSF approximation relevant 
arbitrarily close to the LR, simply by taking v to be sufficiently small. This formal argument allows us to use GSF 
data to explore the immediate vicinity of the LR. 

The structure of this paper is as follows. We start, in Sec. [TTJ by reviewing the formal GSF results relevant to our 
analysis, and then present the new sub-ISCO GSF data. The raw numerical data are given in Appendix [X] for the 
benefit of colleagues interested in reproducing our analysis or studying other applications. In Sec. |III| we use the GSF 
data to construct a global analytic fit for the function a (it), and in particular establish the behavior of this potential 
near the LR. In Sec. IV we similarly construct global analytic models for the 0(v) pieces of £ and J. In Sec.JV] we 
revisit the problem of determining the 0(v) shift in the ISCO frequency, and, using the method proposed in |22j with 
our new, highly accurate a(u) data, add 4 significant digits to the value obtained in previous analyses [TBI [T71 12"2"] . 



Section VI 



turns to discuss the determination of the second 0(v) EOB potential, d(u): by combining the new analytic 
a (it) model with our previously obtained data for p(u), we determine d(u) (numerically) on the domain < 11 < 1. 
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Section VII then focuses on the LR behavior, explaining the physical origin of the observed divergent behavior of 
a(x), and discussing its consequences for the EOB formalism. We summarize our main results in Sec. VIII and discuss 
future directions. 



A. Setup and notation 

Henceforth, we shall use units such that G = c = 1. We will consider a circular-orbit binary of black holes with 
masses m\ < m 2 . Various combinations of these two masses will become relevant in different parts of our analysis: 
we shall use 



M = mi + TO2, 



m 2 



< 1, 



TO1TO2 



(m a 



m 2 ) 



(1) 



to denote, respectively, the total mass, "small" mass ratio, and symmetric mass ratio of the system. This mass 
notation differs from the one used in our previous paper |18j , and is more in line with the notation commonly used in 
EOB and PN work. It reflects the fact that in these formulations (unlike in GSF work) the two masses are treated 
symmetrically. 

We will find it convenient, in different parts of the analysis, to use different measures of the binary separation. 
In the GSF-relevant limit q — > (<^ v — > 0) we will use the standard (areal) radial coordinate r associated with 
the Schwarzschild geometry of the black hole with mass 7712, while in discussing EOB we will mainly use the EOB 
"gravitational potential" (or "inverse radius") u = M/t^ob- A relation between the GSF and EOB descriptions can 
be established using the invariant frequency fl associated with the orbit, or the dimensionless frequency parameter 



(MO) 



2/3 



[(mi 



m 2 



)0] 2 / 3 



(2) 



derived from it. As is well-known, in the GSF limit v — > 0, x becomes equal to u ("Kepler's third law"): x = u + 0(v). 
When discussing the behavior near the (unperturbed) LR, x = | = it, it will be convenient to introduce the (invariant) 
coordinate 



z = 1 — 3a;. 

(The quantity should not be confused with Z\, denoting the redshift of worldline 1. 
summarizes our notation for various mass and radius quantities. 



(3) 

For easy reference, Table [I] 



Binary masses 


Measures of binary separation 


mi 
m 2 

M = mi + m,2 
q = mi/ra2 

— mim 2 
(m 1 +m 2 ) 2 
_ m 1 m 2 


particle mass 
black hole mass 
total mass 
"small" mass ratio 
symmetric mass ratio 
reduced mass 


r or ro 
u = M/r EOB 

n 

x = (Mn) 2/3 

z — 1 — 3x 


Schwarzschild radial coordinate 
EOB "inverse radius" coordinate 
invariant orbital frequency 
dimensionless frequency parameter 
invariant "distance" from light ring 



TABLE I: Various mass and separation quantities appearing in our analysis, summarized here for easy reference. 

II. CONSERVATIVE GSF FOR (STABLE OR UNSTABLE) CIRCULAR ORBITS 

A. Redshift function and regularized self-metric perturbation 

The GSF formulation stems from a perturbative treatment of the binary dynamics. At the limit q — > the object 
with mass mi becomes a "test particle" and its motion is described by some geodesic in a "background" Schwarzschild 
geometry of mass m.2- Finite-mi effects (self- force, including radiation reaction, etc.) are incorporated, in principle, 
order by order in q, working on the fixed background of the large black hole. In this treatment, the small object 
experiences a GSF caused by an interaction with its own gravitational perturbation, and giving rise to an accelerated 
motion with respect to the Schwarzschild background. The GSF accounts for the dissipative decay of bound orbits, 
as well as for conservative (e.g., precessional) effects associated with the finiteness of mi. While the GSF itself is 
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gauge-dependent, knowledge of the GSF (in a particular gauge) together with the metric perturbation due to mi (in 
that same gauge) gives sufficient information for quantifying the gauge-invariant aspects of the dynamics. At the 
foundational level the GSF is now well understood at the first order in q beyond the geodesic approximation [2"5"H2T)] . 
and at this order there is also a well-developed methodology and a toolkit for numerical computations, at least in the 
case of a Schwarzschild background [T^l [3"U] . (The foundations for the second-order GSF have also been laid recently 
|31H33j but this formulation is yet to be implemented numerically.) 

In the problem at hand we ignore the dissipative effect of the GSF, and the orbit is assumed to be precisely circular. 
We shall assume, without loss of generality, that the motion takes place in the equatorial plane 9 = tt/2, where 
hereafter we use standard Schwarzchild coordinates {t, r, 9, ip} defined with respect to the background Schwarzschild 
geometry with metric 3^(7712). Detweiler and Whiting have shown [34] that the GSF-corrected worldline has the 

interpretation of a geodesic in a smooth perturbed spacetime with metric g a p — (m 2 ) + h% p , where h R p (the "R 
field") is a certain [0(g)] smooth perturbation associated with m\. We let uf = {u\,0,0,u±} be the four-velocity 
defined with respect to proper time along this effective geodesic. It is straightforward to show that both u\ and uf are 
invariant under gauge transformations with generators £ Q = O(g) that respect the helical symmetry of the perturbed 
spacetime 35 . The azimuthal frequency (with respect to a coordinate time t belonging to an "asymptotically flat" 
coordinate system), 



n 



dip 



(4) 



is thus also invariant under such gauge transformations. Detweiler |21[ 136) proposed utilizing the functional relation 
u\(n), or, equivalently the "redshift function" 



«i(n) = i/i4(fi), 



(5) 



as a gauge-invariant handle on the conservative effect of the GSF in circular motion. He also discussed the physical 
meaning of z\ as a measure of the (regularized) gravitational redshift between the worldline of m\ and infinity. 

The expressions derived by Detweiler for u\ (fi) [or z-y (fi)] involve the double contraction of h^p with the four- velocity 
u", namely 



h R,G = ,r,g a B 

!h uu — a aB "l "1 • 



(6) 



Here we have introduced the extra label G, for "gauge" (besides the first label R referring to "regularized"), to keep 
track of the coordinate gauge in which one evaluates the metric perturbation. This is important for the following 
reason. The prescription in Ref. |21j assumes that the metric perturbation is given in a gauge which is manifestly 
asymptotically flat (i.e., one in which the unregularized metric perturbation vanishes at infinity). This, however, 
happens not to be the case for the Lorenz gauge that we shall use in our actual GSF calculations. As a consequence, a 
certain gauge correction term will enter our expressions for u\(fl), as we discuss below. We shall use the label G = F 
to refer to a manifestly asymptotically-flat gauge, and the label G = L for the Lorenz gauge. 
Using an asymptotically flat gauge, and the dimensionless frequency parameter 



Ref. [5T] obtained the simple relation 



1 



l + -h 



R,F 



In terms of the redshift variable ([5| this reads 



Zl 



,R,F 



0(g 2 



0(g 2 ) 



(7) 



(8) 



(9) 



The GSF-adapted frequency parameter y [Eq. ([7])] is related to the more symmetric (EOB-adapted) frequency pa- 
rameter x [Eq. pj] through 



V 
x 



m 2 
mi + m 2 



2/3 



1 aW = 1 -h + ^ 



(10) 
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so that 



v/T^y = VT^te + q 7 == + 0(q 2 ) 



Substituting in Eq. ^ then yields, through O(q), 

z\{x) = \J\ — 3x 



(ii) 



(12) 



The form of the last relation is invariant under gauge transformations within the class of asymptotically flat (and 
helically symmetric) gauges. However, our GSF calculations will be carried out in the Lorenz gauge, in which the 
metric perturbation h^ v turns out not to decay at infinity (its monopolar piece tends to a constant value there [37]). 
We need to have at hand the link between the normal "asymptotically flat" h R ^ F and its Lorenz-gauge counterpart 
^uu- The issue was discussed in Refs. [IE1I2S], and we recall here the end result. 

A simple gauge transformation away from Lorenz into a corresponding asymptotically flat gauge is obtained by 
rescaling the Lorenz-gauge time coordinate ti, using 



t F = (1 + a) t L , 



(13) 



with 



VI -3a; 

This defines an F-gauge with metric perturbation given [through 0{q)\ by 

2m 2 



ftoo(r) = *oo(r) + 2a 1 



(14) 



(15) 



with h F p 



h L 



for all other components. Since the gauge transformation relating h F p to h^p is regular, the 



corresponding regularized fields h^'p F and h R p are related to one another in just the same way (this comes from a 
general result derived in [35] )• Evaluating on the mi worldline and contracting twice with the four- velocity, one then 
finds 



,R,F 



R,L 



12a 1- 



2m 2 



I \2 



K) 



which reads explicitly [using (u\) 2 = (1 — 3.x) 1 + 0(q)] 



R,L 



h H ^ = h 



2q 



x(l - 2x) 
(1 - 3a;) 3 / 2 



Inserting Eq. (17) into Eq. (12) finally leads to an expression for Z\{x) in terms of h 



R,L. 
uu ■ 



z\(x) = \1 — 3a; 



1 1 h R.L n x ( l - 

2 uu q (1 - 3a;) 3 / 2 



1 -3a; 



(16) 



(17) 



(18) 



In the above expressions we have not specified the argument in terms of which should be expressed, or — more 
precisely — the specific orbit along which /i^ G should be evaluated. The explicit GSF computations presented below 
actually give h^ L along an unperturbed, geodesic orbit, parametrized by the unperturbed Schwarzschild-radius variable 
m2/V. However, to leading order in q we have TO2A = V + 0(q) = x + 0(q), so that in Eqs. (12) and (18) we can 
simply replace h^^(m2/r) — > h^^(x) [as, of course, h R £* itself is already 0(q) and in our analysis we ignore terms 
of 0(q 2 ) or higher]. 

Finally, we note that h R ^ describes a purely conservative effect of the GSF (even though in practice we shall 
extract h R £* from the retarded metric perturbation). To see this, it is enough to recall Eq. (fSJ) , which relates h R £* to 
the time- symmetric function u\(iY). The property that h R ^ encodes a purely conservative piece of the GSF is special 
to circular orbits, and it does not carry over to (e.g.) eccentric orbits; cf. [4"T] . 
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B. Mode-sum computation of h„. 



Our method follows closely the standard strategy of mode-sum regularization [TJ |39l [40] . As, in this section, we 
work only with the Lorenz-gauge perturbation we shall drop, for concision, the extra label L on h a p. We begin by 
writing = h a p — h^p, where h a p is the full (retarded) Lorenz-gauge metric perturbation associated with the 
mass mi, and h^p is the locally defined Detweilcr- Whiting Singular field ("S field") [34]. Both h a p and h^p diverge 
at the particle, but their difference h^p is perfectly smooth. We formally construct the fields h uu = h a pufu^ and 
h^ u = h^pufuf, where uf is any smooth extension of the four- velocity uf off the particle's worldline (so that uf = u" 

on the worldline itself). We then consider the formal decomposition of the fields h uu (t,r,8,ip) and h^ u (t,r,9,ip) in 
scalar spherical harmonics Y lm (6,(p), defined as usual on the spherically symmetric Schwarzschild background, and 
we let h l uu {r) and h^£(r) denote the individual Z-mode contributions to the respective fields, summed over m for fixed 
Z, and evaluated at the particle [i.e., in the limit r — > r par ti c i e (i)]. As shown in Appendix D of Ref. |41j . the particle 
limit in the above procedure is well defined, and the resulting values h l uu (r) and h^(r) are finite and do not depend 
on the direction (upwards or downwards) from which the limit r — > ^particle (i) is taken. We thus have 

oo 

hL(r)=T,( h uu(r)-&)), (19) 

where it should be noted that while each of the individual Z-mode sums J^i h U u an< ^ ^«« would be divergent, the 
mode sum of the difference ^2i(h l uu — h^) converges exponentially fast (because the difference h a p — h^p is a smooth 
function). We also note that the individual Z-mode contributions h l uu and depend on the off- worldline extension 



chosen for uf, while the sum over modes in Eq. (191 is, of course, extension-independent. 

The formulation of the Z-mode method proceeds by obtaining an analytic description of the large-Z behavior of h^" ■ 
Ref. [5T] (see also Ref. [UJ) obtained the asymptotic form (as I 3> 1) 

h s u i(r) = D (r) + O(r 2 ), (20) 

where D$(r) is an /-independent parameter depending only of the orbital radius r: 

= 4miZ(r) EUipK(w(r))j (21) 
irr 

with 



«(r) = -^-, (22) 
V r — 2rri2 r — Im-i 

and with EllipK(w) = Jq^ 2 (1 — w sin 2 x)~ 1 ^ 2 dx denoting the complete elliptic integral of the first kind. [The large-Z 
behavior of was analyzed in [UJ for generic (stable) eccentric orbits, and we specialize the expressions obtained 
there to circular orbits; the validity of these analytic results for r < 6m_2 will be discussed below.] It was found 
[4T] that the asymptotic value Dq does not depend on the u a -extension involved in the definition of the modes . 
Furthermore it was found that (for any such extension) 

oo 

J2( h uu(r)-Do(r))=0. (23) 

1=0 



This allows us to write Eq. ( 19 ) in the form 



h*u(r) = Y,(h l uu (r)-D (r)), (24) 



which is an operational mode-sum formula for h^ u , describing the correct mode- by- mode regularization of the fields 
h uu . The latter are to be provided as input, typically in the form of numerical solutions to the mode-decomposed 
Lorenz-gauge metric perturbation equations with suitable "retarded" boundary conditions (details of our particular 
numerical implementation arc provided below). 



Since the mode sum in Eq. (19) converges faster than any power of l/l, it follows that the retarded modes to o m ust 



have the asymptotic form h uu (T^> 1) = Dq + 0(1 ). Thus, in general, we expect the partial mode sum in Eq. (24) to 



converge with a slow power law ~ / _1 (and this was indeed confirmed numerically in [SI]). This is problematic from 
the practical point of view, and restricts the accuracy within which h^ u can be computed. As emphasized notably 
in Ref. [ST] the problem can be mitigated by including higher-order terms in the large-/ expansion of h%£- Recently, 
Hcffcrnan et al. [U] were able to obtain analytic expressions for a couple of these: 



Kt(r)=D (r) + 



D 2 (r) , D 4 (r) 



+ o(r 6 ), 



(25) 



where [HHH] L 2 = Q ~ W ~ 
coefficients Z?2,4 are given by 



2)5 ^4 



(J- 



and where the /-independent (but r-dependent) 



D 2 (r) = 



mi 



2Trr 2 Z(r) 



7r 2 — 61m 2 r + 96m 



2m 2 



2 1 EllipK(w(r)) 



(7r - 33m 2 )EllipE(w(r)) 



(26) 



D 4 (r) 



3mi 

V 

1607rr 3 Z(r)(r — 3m 2 ) 2 

30r 5 - 2683m 2 r 4 + 30741mir 3 - 131855m?r 2 + 241905m2r - 160530m§\ „„. T „ , „ 
r 2m2 2 » j E lbpK W r)) 

15r 5 - 1469m 2 r 4 + 13990mir 3 - 56858m3r 2 + 106395m 4 r - 71385m! , 
r-frn, " " ' ™^(Mr)) 



(27) 



with EllipE(ui) = Jq^ 2 (1 — wsin 2 xY^dx denoting the complete elliptic integral of the second kind. [The expressions 
in [42] were derived for generic (stable) bound geodesies, and we specialize them here to circular orbits.] It is important 
to note that the values of the sub-leading parameters D 2 and D4, unlike that of Do, do depend on the off-worldline 
extension of uf. The above values correspond to the particular extension = u i ( m Schwarzschild coordinates), 
i.e., an extension in which the contravariant Schwarzschild components of the field uf are taken to have the constant 
values u" everywhere. This is a practically useful extension and we shall refer to it as the "constant" extension. 

The /-dependent factors in Eq. ( 25 1 have the important property (first exploited in Ref. [13] in the context of the 
scalar-field self-force) 



1=0 A 



00 



(28) 



which allows us to recast the mode-sum formula ( 24 ) in the more useful form 



1=0 



£> 2 (r) D 4 (r) 

Lr> La 



(29) 



Once again, since the sum in Eq. (19) converges faster than any power of 1//, we have that h l uu and must 



share the same asymptotic power-law expansion (25), with the same coefficients D n (as long as h uu is defined and 
computed using the above "constant" u"-extension). Therefore, we expect the revised mode-sum formula ( [29| to 
converge like ~ Z -5 — significantly faster than the original mode sum (24). This will be confirmed numerically below. 



The fast-converging mode-sum formula ( 29 1 forms the basis for our numerical implementation in this work. 



C. Behavior of the mode sum near the light ring 

The results presented in the previous subsection were derived in [HI H2] for stable geodesic orbits. However, all 



of these results, and in particular the form of the mode-sum formula (29) and the values of the parameters D n , are 
equally applicable for circular (timelike) geodesies below the ISCO. Subtleties begin to manifest themselves only when 
the orbit is sufficiently close to the LR at r = 3m 2 . There, the orbit becomes asymptotically null and beaming-type 
effects distort the usual /-mode distribution, potentially enhancing the relative contribution of higher multipoles [see, 
e.g., Davis et al. [44], but note that their analysis concerns the distribution at infinity of £ensoria/-harmonic modes, 
while ours involves sca/ar-harmonic modes near the mi worldline of the particular (extension-dependent) contraction 



9 



That the Z-mode behavior becomes subtle near the LR is evident from the asymptotic form of the parameters D n . 
Defining z = 1 — "imijr we find 

D (z« 1) = - 2qzl/2l f z/16) +0( Z ^ln Z ,z^ (30) 

V37T 

D 2 ( Z « i) = ^- 1/2 [i + m(3z/i6)] + 0(zl/2 lnz , (31 ) 

3V37T 

D 4 (z « 1) = 4q ^_ (z- 7 ' 2 + 13z- 5 / 2 ) + 0(z- 3 / 2 luz, z^' 2 ). (32) 

This suggests that successive terms in the /-mode series become increasingly more singular in 1/z. Even though it is 
not possible to predict the leading-order singular behavior of an arbitrary term Dm based only on the known terms 
-Do, 2, 4 (and this behavior may anyway depend on the extension), it is clear that the limits I — > oo and z — > are not 



interchangeable, and that the mode-sum series (29) becomes ill-convergent near the LR. For any given < z -C 1, 
we expect the series to start showing the standard power-law convergence only for I > Z(z), where l(z) is some 
monotonically increasing function of 1/z, with l(z) — > +oo for z — > + . Our numerical experiments confirm these 
expectations and suggest l(z) oc 1/z — see Fig.[l] 

The evident broadening of the /-mode spectrum near the LR is problematic from the practical point of view: at a 
given z, one must compute at least l(z) modes in order to reach the power-law "tail" regime where the series begins 
to converge, and this quickly becomes computationally prohibitive as z gets smaller. Assuming the empirical scaling 
Z oc 1/z holds, we find that at least ~ 1/z modes must be calculated. Current codes cannot in practice compute more 
than a hundred or so modes, which, a priori, restricts the reach of our analysis to z > 0.01. 

We should comment, in passing, about a more fundamental issue. Strictly speaking, for any (small) nonzero value 
of the mass ratio q, the GSF approximation itself ceases to be meaningful sufficiently close to the LR at z — > 0. This is 
because, for a given q, there are sufficiently small values of z for which the mass-energy of the small particle becomes 
comparable to that of the large black hole, at which point perturbation theory clearly breaks down and the notion of 
GSF is no longer useful. However, reversing the argument, it is also true that we can make the GSF approximation 
valid arbitrarily close to the LR simply by taking q to be sufficiently small. Thus, GSF calculations (and ours in 
particular) can be used to explore the geometry arbitrarily close to the LR. 



D. Raw numerical data for h^ L 

We computed h^ L for a dense sample of orbital radii in the range 3m2 < r < 150to2 using two independent nu- 
merical codes. The first code, presented in Ref. |19j . is based on a direct time-domain integration (in 1+1 dimensions) 
of the metric perturbation equations in the Lorenz gauge. The second code employs a newer algorithm based on a 
frequency-domain treatment of the Lorenz-gauge perturbation equations |45j . E ach code takes as input the orbital 



radius r, and returns the value h^ L (r) computed via the mode-sum formula (29). We typically compute numerically 
the contributions from the modes < Z < 80, confirm the expected Z -6 falloff of the regularized modes (see Fig. [l]), 
and analytically fit a power-law tail to account for the remaining modes 81 < Z < oo. Note that the observed Z~ 6 
behavior comes as a result of a delicate cancellation of as many as six terms in the 1 jl expansion of the unregularized 
modes h uu (r) [i.e., the terms of 0(1°) through 0(l~ 5 )]. It thus provides an excellent cross-validation test for both our 
numerical computation and the analytical parameter values derived in Ref. |42j . 

The new, frequency-domain, Lorenz-gauge algorithm offers significant computational savings as it only involves 
solution of ordinary differential equations, and since, in our circular-orbit case, the spectrum of the perturbation 
fields is trivial (it contains only one frequency for each azimuthal m-mode). This is a crucial improvement, because 
self-force calculations in the time-domain are extremely computationally intensive. The new, frequency-domain code 
allows us to obtain very accurate results at relatively small computational cost. Nonetheless, we have also used our 
time-domain code to check (with lower accuracy) many of our data points. 

Our raw numerical data for h^ L (x), which form the basis for our analysis, are presented in Appendix [a] The data 
for x > 1/5 are new, while our data for x < 1/5 are much improved in accuracy, and more finely sampled, compared 
to previous results [21] [35] . For most data points the fractional accuracy of our data is around ~ 10~ 10 , decreasing 
to ~ 10 -9 at large r and to ~ 10~ 3 very near the LR. (The results of Ref. [24], obtained by a frequency-domain, 
Regge- Wheeler-gauge method, are more accurate than ours, but the data shown in that paper are restricted to the 
weak-field domain 1/500 < x < 1/200.) 
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FIG. 1: Broadening of the Z-mode spectrum near the lig ht ri ng (LR), z = 1 — 3rri2/r — 0. We plot here the (absolute values 
of the) regularized modes h l uu - D - fa - [see Eq. 1 29 1] for < I < 80, for a range of radii on and below the ISCO. The 



dashed lines are arbitrary oc I references. Away from the LR, the regularized modes are expected to fall off at large I with 
an ~ l~ 6 tail, as is clearly manifest in the case z = 1/2 (the ISCO, lower curve). As the radius gets closer to the LR, the 
onset of the l~ 6 tail shifts to larger I- values, with the standard tail not developing until around I ~ 1/z (the regularized mode 
contributions turn from negative to positive around that value of I). In the near-LR case z = 1/200 (upper curve) no transition 
to power law is evident below I = 80. How these data are obtained is described in Sec. |II D| 

Our h^ L data are plotted in Fig. [2] as a function of x. The inset, showing h^ L as a function of z = 1 — 3a; on a 
log-log scale, suggests the near-LR power-law behavior 



■R,L 



(z~ 3/2 asz->0, withC^l- 



(33) 



We will return to discuss the LR behavior in detail in Sec. IVIII 



III. DETERMINING THE EOB POTENTIAL a(x) BELOW THE ISCO 

A. a(x) from h^' F or h^:, L 



Barausse, Buonanno and Le Tiec [53] (using the previous results of Le Tiec et al. 
link between the 0{v) piece zsf{x) of the function z±(x), defined through 



Z\ (x) — vl — 3x + v zsf(x) + 0{v 2 ) , 
and the 0(v) piece a(u) of the EOB function A(u; v), defined through 

A(u; v) = l- 2u + va{u) + 0(v 2 ) , 

namely 

1 - Ax 



i(x) = y/l — 3x zsf(x) — x ( 1 H — I 



]) have derived a simple 
(34) 

(35) 
(36) 
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FIG. 2: Our raw numerical data for the Lorenz-gauge quantity h^ L (x) (blue data points; the solid line is an interpolation). 
The numerical values are tabulated, with error bars, in Appendix [X] Note in the main plot the orbital radius increases to the 
left; the locations of the geodesic ISCO [x = 1/6) and LR [x = 1/3) are marked with vertical dashed lines. The inset shows the 
same data (in absolute value) plotted against z — 1 — 3a; on a double logarithmic scale (note here the orbital radius increases 
to the right, and the LR limit is z — > asymptotically far to the left). The dashed (magenta) line is a simple power-law model 

f,R,L 9 7 -3/2 

'''uu 2 



Let us clarify that in Eq. (36) we regard a(-) and Zsf(') as functions , with x merely denoting a "dummy" argument. 
We could as well have written (36) using the natural notation for the EOB argument of the function a(-), namely 



i(u) = y/l - 3lt z SF (u) - u (l + ~^j=^) 



(37) 



In the following, we shall freely alternate between using x or u as independent variables for the function a(-). Note 
that the two physical variables x and u satisfy the functional relation x(u) — u + 0(v) (see below), so that Eqs. (36) 



and would anyway have the same content (at leading order in v) even if we interpret the arguments x and u as 
physical variables. 



Putting together the definition (34) of zsf(x) and the previous links (12) and (18) between the redshift function 
Z\(x) and the metric perturbations, we have the following relations between zsf(x) and the two types of metric 
perturbations (Flat or Lorenz): 



zsf{x) = y/l - 3z 



2 uu 



1 - 3x 



(38) 



ZSF 



(x) — Vl — 3x 



1 IR,L _ x C 1 ~ 2x ) , x 

2 uu (1 - 3x) 3 / 2 1 - 3x 



Here, the tilde over /i^ indicates that one has factored out the mass ratio q = v + 0(v 2 ): 



~ h R,G = -l h R,G 



Finally, inserting Eq. (38) or Eq. (39) into Eq. (36) yields, respectively, 

1 _ , 7 RF 1 - 4:X 



a(x) 



1 - 3x) h 



VI -3a;' 



(39) 



(40) 



(41) 
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a(x) = -~(l-3x)h%'-2xJT=te. (42) 

Note in passing that some cancellations took place when replacing zsf{x) in terms of either h^ F (x) or h^ L (x). 
In particular, when relating a(x) to the Lorenz-gauge perturbation h^ L , the single "extra" term that remains (the 
one not involving J&P) has the property of tending towards zero at the LR limit (x — > |), while the corresponding 



extra term in Eq. pd|) tends to infinity in that limit. We will come back later to a deeper discussion of the physics 
near the LR. 



The "doubly rescaled" function &e(x) 



Using Eq. (42 1 and our numerical results for h^ L (x), one obtains a dense sample of numerical values for a(x) over 
the entire range < x < |. Our goal now is to obtain a global analytic fit formula that faithfully represents this 
relation. 

The function a(x) itself varies rapidly at both ends of the above domain (just like h^ L in Fig. [ij: as we discuss 
below, a(x) vanishes fast at x — and blows up at x = 1/3. Rather than fitting directly for a(x), it is more 
convenient to fit for a new function, constructed from a(x) by "factoring out" suitable terms representing the leading- 
order behavior at both ends of the domain, so that the resulting "rescaled" function is relatively slowly varying over 
the entire domain. 

Let us first consider the behavior near x = 0. Information from PN theory determines the form of a(x) in this 
weak-field regime. Refs. [U HH [23j EH |46l [47] obtained the expansion 

oo 

a PN (x) = ^2(a n + a l £\nx)x n , (43) 

n—3 

where alf = aq* = 0, and the first few nonzero coefficients that can be determined analytically from available PN 
expressions are 

a 3 = 2, 

94 41tt 2 



In 



3 32 
64 

T' 

7004 



105 



(44) 



Note that the leading-order (Newtonian) behavior is a(x) ~ 2x 3 , and that logarithmic running first appears at 0(x 5 ^ 



A few more, higher-order coefficients in the expansion (43) were obtained numerically in Ref. |23| by fitting to the 
accurate large-radius GSF data of Ref. [2H US] : 

a 5 = +23.50190(5), 
a 6 = -131.72(1), 
a 7 = +118(2), 

4 n = -255.0(5), (45) 



where in each case a parenthetical figure indicates the estimated uncertainty in the last decimal place. Let us introduce 
the notation a(x) to denote the normalization of the function a(x) using the leading-order PN term, i.e., 

2W = ^. ( 46 ) 

so that d(0) = 1. 

Consider next the behavior near the LR, xlr = |. This behavior has not been studied so far, neither numerically, 
nor analytically. Note, however, that Ref. [55] has remarked that the extrapolation beyond x = | of a 5-parameter 
fit to 55 data points (ranging between x = and x = h) for zsf(x) indicated the possible presence of a simple pole 
Zsf{x) ~ (x — a^poic) -1 located near the LR. (Their fit yields £ po i c ~ 0.335967, which is slightly beyond the LR.) In 
this work, we study the behavior near x = | both numerically and analytically. We have already mentioned that 
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our data suggest the scaling relation (|33| (which corresponds to a simple pole in zgp (x) located exactly at the LR) 



Combined with Eq. ( 42 ) , this suggests the leading order divergent behavior 



a{x^\)~\{l-Zx)-V\ (47) 



where, recall, the "fudge" factor £ is ~ 1. We will discuss the analytical origin of this asymptotic behavior in Sec. |VII| 
below. 

Equation (47 1 suggests that it would be convenient to further normalize the function a(x) by a factor (1 — ix)^ 1 / 2 . 



However, as will be further discussed below, there is a more physically motivated normalization: we recall that 
the conserved specific energy associated with mi as it moves along a circular geodesic orbit in the Schwarzschild 
background of ttt-2 is given by 

£ <*> = TO < 48 > 

which has the same type of divergent behavior as a(x) for x — > | (but is regular elsewhere). We hence choose to 
use E for our second normalization. Let us introduce a notation whereby a sub-index E denotes normalization with 
respect to E(x); in particular, 

a E (x) = 49 
E(x) 

Note a B (l/3) ~ f C 

Let us finally introduce the "doubly-rescaled" function 

^ s Jwr < 50 > 

which attains finite, nonzero values at both ends of the domain < x < |, namely o,e{0) — 1 and g,£(1/3) ~ ^C- 
Our numerical dataset for cle{x) is plotted in Fig. [3] Evidently, the function £le(x) is monotonically increasing and 
convex over < x < |. 



C. Accurate global analytic model 

We have explored a large set of global analytic models for the function cle(x). In each case we used a least-squares 
approach, i.e. given the data (sampled at a discrete set xi,X2,--- of x values), together with an estimate of the 
corresponding standard data errors <7d a t a (xj), we minimized (using Mathematical function NonlinearModelFit [ ] ) 
the standard x 2 statistic 

9 . , v-^ /data(xi) — model(.Ei, parameters) \ 2 

X £ parameters = > — (51) 

^ \ cr data (a:j) J 

over the model parameters to determine the best-fit parameters. In addition, we evaluated the faithfulness of the fit 
by recording the minimum value of \ 2 '- Xmm = X 2 (best fit parameters) . If numerical errors were normally distributed 
without a systematic bias (which we assume here) , and if the model were "true" , then we would expect Xmin to equal 
approximately the number of degrees of freedom (DoF) in our model. As a second measure of the fit quality we also 
considered the norm \ \a^ — a^f ta || , i.e., the maximal absolute difference between the best-fit model and the data 
over all data points. We are using here a^(a;) [rather than a{x) or cle(x)} because this is the relevant quantity entering 
the EOB expressions, as we discuss in Sec. |VII| below. 

We report here some of our results, and present two selected models: a 16-parameter high-accuracy model with 
% 2 /DoF of order unity; and (in the following subsection) a simpler, 8-parameter model, which is less accurate (has a 
very large x 2 value) but has a sufficiently small norm | [a^f — a|f ta | | to be useful in some foreseeable applications. 

Let us focus, for ease of presentation, on the following restricted class of analytic models, which employ a Pade-like 
approximant for ag: 

1 + ELi id + c[° S Inx) x i + x 3 z\n \z\ (c zQ + c l °§ In \z\ + c zl z) 

4* = — ^ \,^ q , ,■ L - (52) 

1 + ELi djxJ 
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FIG. 3: Numerical data for the doubly-rescaled function &e{x) [see Eq. ( 50 1] . The solid line is a cubic interpolation of the 
numerical data points (beads). The inset shows, on a semi-logarithmic scale, the relative numerical error in the &e data, 
computed based on the estimated errors tabulated in AppendixlAl Note that the relative error is between 10 -8 and 10 -10 over 
most of the domain, and it never exceeds 10~ 5 (except at a single point, closest to the LR, where it is ~ 0.1%). 



Here p > 3 and q>l are constant integers (see below), and {ci>2, c i>& di> c z0> c z0-> c zi} are the model parameters to 
be fitted. The first few c parameters are constrained so as to reproduce all analytically available PN information: 



97 

o, 

32 



4l7T 2 

~64~ 



3166 32 



105 



di. 



(53) 



We do not constrain the remaining parameters to agree with the additional PN information available through numerical 
fit [Eq. (45 1], but rather allow our model to "re-fit" some of these high-order PN terms. We find, in general, that this 
leads to improved global fits. 

Our model family a^{x) is designed (heuristically) to capture all global features of cle(x) from x = down through 
the LR and (potentially) beyond. We use a Pade-type expansion in x (with logarithmic running terms), augmented 
with z-dependent terms which are aimed at capturing the behavior near the LR. The latter terms are multiplied by 
x 3 to suppress their support in the weak field, where the known PN behavior should apply (we have tried various 
powers of x and found that x 3 generally works best). 

The form of the ^-dependent terms in (52 1 is motivated as follows. We have initially experimented with simple 



polynomials in z (without logarithmic terms), but found that these always yielded best-fit models that possessed poles 
(singularities) immediately behind the LR (i.e., just above x = |). This suggested to us that the true function cle(x) 
has a remaining non-smoothness at x = |, and the form of the function suggested a weakly divergent derivative. In 
our model family (52 1 we have attempted to represent this type of non-smoothness with a term of the form ~ zLn|z|, 
which indeed seemed to have the effect of removing the undesired pole. To allow more freedom in fitting the correct 
LR behavior we have added a few higher-order terms in z and In \ z\. We experimented with a large variety of such 



higher-order term combinations, and found that the form shown in (52) worked well (while minimizing the number 
of extra model parameters). 

Each member of our model family a^(x;p, q) has 2p + q — 1 fitting parameters. In Table [TT| we show fitting results 
for a variety of p, q values (and also for models in which we remove some of the In x terms) . For each fitting model we 
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model 

# 


Fit model [Eq. (52 1] 


# model 
parameters 


X 2 /DoF 


II of -4? ta || 

II B K II CO 


pole? 


c 


P 


q 


params. set to zero 


1 

2 


4 


4 


C4° S 


10 


4.45 x 10 6 


1.62 x 10" 2 


— 


0.991785 




11 


6.71 x 10 4 


4.01 x 10~ 3 


— 


1.00984 


3 
4 


4 


5 


4 


11 


3.81 x 10 4 


2.65 x 10 -4 


— 


1.00791 




12 


2.81 x 10 3 


1.46 x 10" 4 




1.00192 


5 
6 
7 


5 


5 


c k)g c log 


12 


5.55 x 10 3 


4.53 x 10" 5 


— 


1.00408 


c' og 




13 


2.79 x 10 3 


1.37 x 10" 4 


— 


1.00214 




14 


1.52 x 10 3 


1.48 x 10~ 4 


x « 0.45 


0.999680 


8 
9 


5 


6 


Jog 
5 


14 


1.83 x 10 3 


1.41 x 10" 4 


— 


1.00137 




15 


1.03 x 10 3 


1.98 x 10" 4 


x ~ 0.35 


0.988620 


10 
11 
12 
13 


6 


7 


log log log 
4 ' o ' fa 


15 


19.2 


1.12 x 10" 6 


— 


1.00750 


log log 
5 ' 6 


16 


9.97 


1.08 x 10" 5 


x w 0.42 


1.00536 


log 
6 


17 


3.37 


2.52 x 10~ 6 


x w 0.375 


1.00525 




18 


4.94 


2.01 x 10~ 5 


— 


1.00907 


14 

15 
16 
17 
18 


7 


7 


log log log log 
4 ' o ' fa ' * 


16 


4.77 


1.97 x 10 5 


— 


1.00899 


log log log 
a ' fa ' 7 


17 


3.08 


3.81 x 10" 6 


x w 0.36 


1.00453 


log log 
6 ' 7 


18 


3.03 


5.81 x 10" 6 


x w 0.35 


1.00345 


c' og 


19 


2.87 


4.28 x 10~ 6 


x w 0.575 


1.00968 




20 


2.93 


2.62 x 10" 6 


x w 0.58 


1.00918 


19 
20 
21 
22 
23 


7 


8 


log log log log 
c 4 i L 5 i 6 i c 7 


17 


4.79 


1.79 x 10" 5 


— 


1.00882 


log log log 
L 5 ' 6 i L 7 


18 


3.04 


5.71 x 10" 6 


x « 0.35 


1.00350 


log log 
L 6 ' L 7 


19 


3.01 


5.46 x 10" 6 


x w 0.35 


1.00373 


Jog 

c 7 


20 


2.87 


2.51 x 10 -6 


x w 0.53 


1.00775 




21 


2.87 


2.16 x 10" 6 


x « 0.5 


1.00732 


24 


4 


4 


log log 
°4 ' c 'z0 ' 02:1 


8 


1.08 x 10 7 


1.21 x 10" 5 


x ~ 0.7 


1.00554 



TABLE II: Model fitting for the doubly-rescaled function &e(x). Each row describes best-fit results for the model af (x) given 
in Eq. (52 1, with particular values of p and q; in some of the models we have eliminated some of the fitting par ame ters, as 
indicated in the fourth column. [In the last row we show best-fit results for the model a f E t ' s ' mp (x) given in Eq. | |54| , to be 
discussed in Sub. HID below]. The fifth column shows the total number of fitting parameters for each model, and the sixth 
columns displays the value of x 2 /DoF for the best fit parameters. In the seventh column we show the maximal difference 
between the model (with the best fit parameters) and the data, for the physically relevant quantity o,e{x) = a{x) / E(x). In 
the penultimate column we indicate the location of the first pole of af (a:) [w hich is the same as for af (a;) or a Rt (x)}. The 
last column presents the value of the fudge factor £ = |a.e(l/3) [see Eq. ( 47 1] , as predicted by each of the best-fit models. 
Highlighted in boldface are values for the two selected models (#14 and #2^ 
Tables HTT1 and ITVl below. 



whose parameters are given, respectively, in 



compute the x 2 statistic using as weights the estimated numerical errors from Tables VIII and IX For each best-fit 
model we also display in Table [n] the value of the norm 1 1 af — a^f ta 1 1 ^ . Some of the models presented in Table [TT| 
have remaining poles on x > 1/3, and we indicate in the table the location of the first pole below the LR if any occurs. 
Finally, the table shows the predicted value of the fudge factor £ for each of the models. 

The data in Table |Tl| suggest that, at least within the model family (52), one cannot obtain a good fit for cle{x) 
with just a handful of model parameters. At least 14 parameters are needed to achieve x 2 /DoF ~ 1000 and at least 
16 for x 2 /DoF < 10. However, the value of x 2 /DoF becomes saturated at around 3 or 4 for > 16 parameters, and 
does not decrease much further upon adding extra parameters (this may indicate that our quoted numerical errors 
<7i are slightly too optimistic, consistent with our estimated factor ~ 2 uncertainty in the values of the quoted a^; see 
Appendix |A| . We have experimented with several other model families but were not able to achieve x 2 /DoF values 
of order unity with less than 16 parameters. 

We choose to present here the accurate 16-parameter model highlighted in Table [IT] (model #14), which has 
X 2 /DoF=4.77 and ||af — a £; ata || 00 = 1-97 x 10 -5 . Some of the other models in the table have slightly lower val- 
ues of x 2 /DoF but in all such cases the models are more complicated (have more parameters) and, more importantly, 
they present undesired poles below the LR. The best-fit parameters for our selected model are shown in Table |III| 
Note that we are giving the parameter values accurate to many decimal places: this accuracy is necessary for our 
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i 


Ci 


Jog 


di 


1 


+7.48610059021 





-2.357850757006 


2 


-8.81722069138 


, 32 
' 5 


-1.889967139293 


3 


-227.6806641934 


-45.2426257972 


-109.86788081837 


4 


-1336.5402672986 





+535.45853874191 


5 


+8044.588011262 





-53.572041734 


6 


-5643.745303388 





-3030.7781195456 


7 


-7744.83943928 





+4106.962599268 



CzO = 


-32.8395937428 


J°g _ 

C z0 — 


-4.34430971904 


C Z 1 = 


-365.569972774 



TABLE III: Parameter values for the 16-parameter model highlighted in Table \H\ The model belongs to the family (521, with 
ci and c±2 3 constrained in accordance with (531 so as to impose all PN information available analytically. This 16-parameter 
model has x 2 /DoF ~ 4.77 and it reproduces all of our numerical data points for cle(x) to within an absolute difference of 
1.97 x 10 -5 . (The difference for most data points is actually much smaller — see Fig. [4]) All decimal figures are significant, in 
the sense that removal of any figure would lead to X* /DoF > 4.77. 



10" 



10" 



10" 



io" 10 . 



I a.E (model) — cle (datd)\ 



10" 



10" 




0.05 0.1 0.15 0.2 0.25 0.3 0.35 

x 




0.35 



FIG. 4: Faithfulness of the analytic best-fit model (52 1, with parameters as given in Table III The left panel shows, on a 
semi-logarithmic scale, the magnitude of the absolute difference between the model and the data; we use here the variable 
o,e(x) [rather than &e{x)], which is the relevant one entering the EOB potential. The right panel shows (now on a linear scale) 
that same difference divided by the estimated numerical error for each data point. For most data points the model reproduces 
the data down to the level of our numerical noise. 



model to reproduce the data down to the level of the numerical noise, which is generally as small as 10 -9 -10 -10 
(fractionally). We have checked that all decimal figures shown are significant, in the sense that omitting any of the 
figures would result in x 2 /DoF values larger than 4.77. 

Figure [4] shows a performance diagnostic for our selected model. From the left panel it is evident that our model 
for <1e(x) reproduces most data points to within mere differences of 10 -10 -10 -12 . Larger differences appear only at 
x > 0.3, in which domain our data is less accurate. The right panel compares the difference between the model and 
the data with the estimated numerical error in the data points. We observe that most data points are reproduced by 
the model at the level of the numerical noise, as desired. 

In Fig. [5] we plot our selected model for cie{x) over the entire domain < x < 1, showing how it extends beyond 
the LR. We observe that the function ole{x) peaks closely below the LR, then drops and changes its sign around 
x ~ I (the location of the event horizon in the background geometry). To assess the robustness of these features 
we have also plotted in Fig. [5] the global extensions of a handful of other models — all models from Table |III| with 
\\a^ — a^f ta 1 1 < 10~ 4 admitting a smooth behavior with no poles below the LR. Remarkably, the above basic 
features seem to be preserved: a peak right below the LR, followed by a change of sign. The function cle{x) generally 
turns negative at x values in the range ~ 0.5-0.6 (more towards 0.5 for the more accurate models). Note that the 
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model # 4 
model # 5 
model # 6 
model # 8 
model #10 
model #13 
model # 14 
model #19 




FIG. 5: Extension of our analytic cle(x) models below the LR. The thick (blue) curve shows the behavior of our selected model 
(#14 in Table [TT] and Table [TTT| ) over the entire domain < x < 1. Other curves, labelled by model numbers from Table |Il) 
show the behavior of other models for comparison. Shown, from top to bottom at x = 0.8, are models number 5, 8, 6, 4, 13, 
14, 19 and 10. 



function 0^(1) [as well as a(x) itself] turns negative at precisely the same location as Note also that oe(x) 

vanishes at the proximity of x = | despite having E{x) (which vanishes at x = |) being factored out in its definition 
[recall cie(%) = a{x) / E(x)]. This may (heuristically) point to a rather rapid vanishing of a(x) at the horizon. Whether 
or not the above features are indeed robust remains to be verified. 



D. Simpler global analytic model 



In the above "high fidelity" 16-parameter model, the maximum norm of differences determined 
by the data point nearest the LR, which is our least accurate point. As can be seen from Fig. [4] (left panel), removing 
just a few near-LR points from our sample would result in a norm of a mere ~ 10~ 10 . This high standard of accuracy 
may not be necessary for some applications. Indeed, in designing analytical waveform templates for comparable-mass 
binaries one usually has no reason to require the model to be locally accurate at that level. It is therefore both 
convenient and useful to have at hand a simpler, less unwieldy model, which reproduces all numerical data points to 
within a prescribed accuracy (say, a fiducial < 10 -5 ), but not necessarily to within the (very high) accuracy of our 
numerical computation; in other words, a model in which we relax the requirement % 2 /DoF ~ 1 and replace it with 
an upper bound on the norm I la 



at 



-data I 



,fit 



-.data I 



Through experimentation, we were able to devise an 8-parameter model with a norm as small as 
1.2 x 10 -5 (that is, even slightly smaller than the norm for our 16-parameter model). The model is given by 



a 



fit, simp 

E ' ~ 



1 + c\x + x 2 (c2 + c 2 ° s lnx) + x A (cz + c£ s lnx) + C4X 4 + c z qx 4 z In \z\ 
1 + d\x + d,2X 2 + d^x 3 + d^x 4 



log 



(54) 



where C\ and c 2 °| are again PN-constrained as in Eq. (53), and {C2, C3, C4, c z q, di, d,2, #3, 04} are 8 independent model 
parameters. The best-fit values of these parameters are given in Table |TV} and the residual differences between the 
model and the data are plotted in Fig. [6} We see that the model reproduces the numerical as(x) data to within 10~ 7 
for x < 0.28 and to within 10~ 6 for x < 0.32. 

We note that the above simple model has a very large value of x 2 /DoF (~ 1.08 x 10 7 ). Also, its behavior beyond the 
LR is problematic and rather different from that shown in Fig. [5j the function o>e(x) does not reach a maximum but 
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i 






di 


1 


+2.0154525 





-7.82849889 


2 


-39.57186 


,32 

5 


+20.84938506 


3 


-24.30744 


— 80.254774 


-20.5092515 


4 


+103.93432 





+5.383192 



CzO 



40.22474 



TABLE IV: Parameter values for the "simple" model (54 1. The parameters Ci and c^f 3 are constrained in accordance with 



Eq. ( 53 1 , and the other 8 parameters are found by model fitting. This 8-parameter model reproduces all of our numerical data 
points for oe(:c) to within an absolute difference of 1.2 x 10 -5 . All decimal figures are significant, in the sense that removal of 
any figure would lead to a larger difference. 



\a,E(mode I) — ajs;(datd)\ 
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FIG. 6: Performance of the simpler analytic model (54 1, with the parameters given in Table IV The plot shows, on a semi 



logarithmic scale, the magnitude of the absolute difference between the model and the numerical data for the function o,e(x)- 



instead it grows monotonically with x and eventually blows up at a pole located at x = 0.695694 . . . (i.e., below the 
background horizon). However, we emphasize, the model reproduces all of the numerical data points for the function 
oe(x) to within a maximal absolute difference of a mere 1.2 x 10~ 5 over the entire domain < x < i. 



E. PN limit of the global analytic models 



It may be useful to study here the PN expansion of our global models. Recall that in our treatment we have imposed 
all PN information known analytically [Eq. (44)] but refrained from imposing the numerical values of the higher-order 
PN coefficients a^,ae, . . . [Eq. (45 1] obtained in [23] by fitting to large-radius numerical data from Refs. 24, 46]. We 
may now check how these numerically specified high-order PN coefficients compare with the ones entailed by our 
global analytic models. By considering the PN expansion of our 16-parameter analytic model (#14), we obtain, for 
the two leading coefficients, 



,#14 
,#14 



= 23.47267... [23.50190], 
= -127.154... [-131.72], 



(55) 



where in square brackets we recall the values from [23J [Eq. ( 45 ) above] . Higher-order coefficients agree only in sign 



i* 14 = 701.092. . . [118] and a 7 n ' #14 - 83.0457. . . [-255.0]. Similarly, our simple, 8-parameter model (#24) yields 



#24 
,#24 



24.19028. 
-163.396. 



[23.50190], 
. [-131.72]. 



(56) 



with a-j and of agreeing much less well (and a-j not even agreeing in sign). 
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It should come as no surprise that the values extracted from our global fits differ from those obtained using large-r 
data only, and that the discrepancy increases rapidly with PN order. Our goal here was not to obtain accurate values 
for the coefficients of the asymptotic PN series (as in Refs. [HI SB]) but rather to devise a globally accurate model 
for the function a(x). The latter goal could be achieved more "economically" (i.e., with a simpler analytic model) by 
relaxing (and thus effectively re-fitting) the values of some of the high-order PN coefficients. 

provides a rough idea of the uncertainty within which 



Finally, comparing between the results of Eqs. (55) and (56 



the values of PN parameters can be extracted from any globa fit. We see that, while different global models roughly 
agree on the value of as, they predict rather different values for a$ (and for 07, a' 7 n , . . .). One should keep this 
uncertainty in mind when comparing the "effective values" of the PN coefficients extracted from any single global fit. 

IV. GSF CORRECTION TO THE ENERGY AND ANGULAR MOMENTUM OF CIRCULAR ORBITS 

A. Energy and angular momentum in terms of a(x) and a'(x) 



Reference |H] has used the results of J2U] to derive links between the function zsf{%), Eq. (34), and the GSF 
corrections to the functions e(x) and j(x), where e = {£ — M)/fj, is the binding energy per unit reduced mass, 
and j = J/(Mfjb) — J/[mim2) is the rescaled total angular momentum. Here £ represents the (invariant) total 
gravitational energy of the binary, as it is defined in PN or EOB theories [55]. By inserting the link po) between 
zsf{%) and a{x) into their results one can derive the corresponding links between (e(x),j(x)) and the function a(x). 
However, we find it simpler, and conceptually more transparent, to use (more general) known results from EOB theory 
to directly derive the latter links. Let us start by recalling some basic results from EOB theory of circular orbits and 
its GSF expansion (see [H] for details). 

The total energy £ (including the rest mass contribution) of a (circular) binary system is simply given by the value 
of the EOB Hamiltonian H. Given the EOB main radial potential A{u\ v), it is easy to derive the exact link between 
H and the EOB radial variable u = M/tbob- (This is done, exactly as in the textbook treatments of circular orbits 
around a Schwarzschild black hole, by extremizing an effective potential; see below.) Before describing the result of 
this extremization let us recall the explicit structure of the EOB Hamiltonian: it is given by 



H(u,j,p r ) = Mh(u,j,p r ) , (57) 



with 



and 



h= a/1 + 2i/(iJoff - 1) (58) 



H cS (u,j, Pr ) = ^A{u-v) (l+fu 2 + +Q(u,j,p r ;^y (59) 

where the second EOB potential B(u; v) = <?°f is related to A and the potential D(u; v) mentioned above via 

ABD= 1. 

Along circular orbits p 2 vanishes, as does Q(u, j,p r ;is) when one is working within the standard formulation of 
the EOB Hamiltonian (see below for more details). Circular orbits are obtained by extremizing the simple effective 
potential A{u;v){l + j 2 u 2 ) with respect to u, at fixed j. This leads to the following relation (valid along circular 
orbits) between j 2 and u: 

■>, x A'(u) A'(u) ,„„. 

f{u) = — 2 ) ' =--4/-, 60 

where a prime denotes differentiation with respect to u, and we have introduced the shorthand notation 

A{u) = A(u) + ^uA'(u). (61) 



[For notational brevity we hereafter ignore the ^-dependence of A(u, v).] Note that Eq. (60) yields the simple result 



(62) 
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Inserting the latter result into the effective Hamiltonian, Eq. (59), yields yet another simple result: 



H e g{u) 



A{u) 



A(u) 



(63) 



At this stage Eqs. (60) and (63) give the exact functional relations between u and the energy and angular momentum 
of circular orbits. Note that the specific binding energy e = (H — M) / fi reads, in terms of the above notation, 



H-M 



v v 



1 + 2v(H cS - 1) - 1 



(64) 



In order to obtain the corresponding exact functional relation between the frequency parameter x and e and j 
we need to relate x to u. The latter link simply follows from the general Hamiltonian equation = dH/dJ, which 
explicitly reads (along circular orbits) 



j(u) u 2 A(u) 
h(u)H eS (u) 



Squaring this result, and substituting from Eqs. (60) and (163]), yields the simpler looking relation 



M z n z (u) = - 



1 A'{u) 

2 h?{u) 



(65) 



(66) 



Recalling the definition of x in Eq. we simply have (Mfl) 2 = x 3 so that (66) yields the following exact link 
between u and x: 



x(u) = u 



-\A'{u) 
h 2 {u) 



1/3 



(67) 



Note that, up to this stage, we have made no approximation. In other words, given an explicit expression for the 
main EOB potential A(u;v), one can, by using EOB theory, write the exact functions e(u), j(u) and x(u) along the 
sequence of circular orbits. Note also that we are considering here the full sequence of stable or unstable circular 
orbits. (See [IB], and below, for the exact condition defining the ISCO separating stable orbits from unstable orbits.) 

As, in this work, we are interested in GSF expansions at order 0(v) = O(q), let us now expand the above results in 
powers of v. Using Eq. (35), it is trivial to obtain the functions e{u) and j(u) to order 0(v) in terms of the function 
a(u) [see, e.g., Eq. (4.19jof 16J for the expansion of j(u)]. However, an extra complication comes from the need to 
also expand the function x(u) to order 0(v). This result was first obtained in Ref. |TB] , Eq. (4.21). Here, we are 
mainly interested in the 0{u) inverse relation u(x), which reads [Eq. (4.22) in |16j ] 



u(x) = x 



1 „ s 2 / l-2x 

l+ r a ^ + r\wTTTx- 1 



0{y 2 



(68) 



Finally, inserting (68) into the 0(y) expansions of e{u) and j(u) leads to the following 0(y) expansions of the 
composed functions e[u[x)) and j(u(x)): 



e(u(x)) = e (x) + v [a E (x) a'(x) + /3 E (x) a(x) + Je(x)] 
= e (x) + ve S F(x), 



(69) 



where 



j(u(x)) = j (x) +v[a j (x)a'(x) +(3 j (x)a(x) +jj(x)} (70) 
= jo 0*0 + vjsF{x), 



1 — 2t 

«W = ^j=g - 1, (71) 



Mx) s 7m=wr (72) 
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and where the coefficients cie, /3e, je, ay , j3j and 7^ are given by 

1 x 

Oe(x) 



3 VI -3a;' 

1 l-4x 

2 (l-3x) 3 / 2 ' 



7b W = - e o(z) x 



1 . . 2; 1 — 6x 
2 eo(x)+ 3 (l-3:r)3/' 



a (x) = x 3/2 a E (x) = -i 1 ■ 

•J \/s(l — 3x) 



2x 1 /2(l-3 a; )3/2' 



7j(a?) 



1 



1 — 6a; 



x e (x). 



(73) 
(74) 
(75) 
(76) 
(77) 
(78) 



3x 1 /2(l-3 ;r )3/2 

Note in passing that the relation de[u(x)} = x 3 ^ 2 dj[u(x)] holds exactly (and, in particular, at each order in u), and 
implies many relations between the various coefficient functions ole\ (3e, 7_e, ctj , /3j, jj. The simplest of these relations 
is the link oe(x) = x 3 ^ 2 aj(x) indicated above. 



B. Global fits for the binding energy and angular momentum 



Equations (69) and (70) link the functions esi?(x) and Jsf{x) to the function a(x) and its derivative a'(x). Our 
global analytic model (52) hence translates to global analytic models for e S F{x) and jsF{x)- Since the resulting 



analytic expressions are quite cumbersome we will not present them here but instead content ourselves with a plot of 
the results. 

First, however, it is useful to consider the asymptotic behavior of esF^x) and jsF(x) at the two ends of the domain 
< x < 1. In the weak- field regime x<l, our models, of course, reproduce the known PN behavior: 

e SF (x« 1) = l 7 x 2 + 0(x% (79) 
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Jsf(x«1)= 1 -x^ + 0(x 3 / 2 ). 



Near the LR, we use Eq. (47) in conjunction with (69) and (70) to obtain 

1, 



zsf{x -> 



1-1 C 
27 12 s 



108 



(80) 



(81) 



Jsf{x -> -) 



3\/3 



c »- 



12-s/3 ' 



(82) 



where we also used 1. Hence, both e$F an d jsF are expected to diverge quadratically in 1/z = (1 — 3x) 1 at the 
LR. 

To plot our analytic models for the binding energy and angular momentum it is convenient to introduce the rescaled 
quantities 

e SF = e SF x z 2 x {x 2 /2A)~ l , (83) 

Jsf = Jsf x z 2 x (Vx/6)-\ (84) 

which attain the regular values csf(0) = Jsf(O) = 1 as well as es_F(l/3) « —10 and jsf(1/3) sa — |. The functions 
csf(x) and Jsf{x) are plotted in Fig. [7] Inspecting the plot, we note that both GSF corrections Ssf(x) and jsF(x) 
turn from positive in the weak field to negative in the strong field. The transition occurs at x « 0.0247 (r ~ 40.49M) 
for e SF and at x « 0.0435 (r « 22.99M) for 
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■ (normalized) GSF correction to angular momentum 
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FIG. 7: The GSF corrections to the binding energy e(u(x)) and angular mom entu m j (u(x)). We show here the rescaled 
functions esF(x) (lower curve) and Jsf{x) (upper curve) defined in Eqs. ( 83 1 and (84 1. The curves shown represent the 
analytic functions obtained by inserting our global analytic fit for a(x), Eq. (52 1, into Eqs. (69 1 and ( TO I . 



V. ACCURATE DETERMINATION OF THE Oiv) CORRECTION TO THE ISCO FREQUENCY 



Le Tiec et al. pointed out in Ref. [55] that the link they had established between the functions e(x) and zsf(x) 
provides an efficient method for calculating the GSF- induced [O(^)] shift in the value of the ISCO frequency — 
an important strong-field benchmark. This shift was first calculated in Ref. [T7] (with a crucial gauge correction 
introduced later in |16j ) by analyzing small-eccentricity perturbations of circular orbits. Full details of this calculation 
(and a slightly more accurate result) were presented in [TJ]]. Le Tiec et al. suggested calculating the ISCO shift by 
minimizing the binding energy function e(x). This seems potentially advantageous from the computational point 
of view, because the function zgp(x) [from which e(x) is determined] is derived from GSF computations along the 
sequence of strictly circular orbits. These are substantially simpler and less demanding than GSF computations 
along eccentric orbits, even for the small eccentricities considered in [T7l[T(5]. [However, below we comment that Le 
Tiec et a/.'s method is essentially equivalent (computationally) to the second method used in |17[ If 9j. in which the 
GSF is computed along circular orbits with a certain fictitious source term containing derivatives of the particle's 
energy- momentum.] 

The calculation by Le Tiec et al. in Ref. [55] , based on the circular-orbit GSF data for r > 5m2 available to them at 
the time, produced a value in full agreement with the results of [HI [TTl [19] , and with a slightly improved accuracy — see 
Table |V} This agreement also lent support to the assertion made in Ref. [50] that the link between e(x) and z$f(x) 
is valid not only through 3PN order (as explicitly proven by them) but also in the strong-field regime. 



Source 


Cn 


Barack & Sago 17 ; Damour |16| 
Barack & Sago [19] 
Le Tiec, Barausse & Buonanno [52] 
This work 


1.2513(6) 
1.2512(4) 
1.2510(2) 
1.25101546(5) 



TABLE V: Value of the parameter Cn describing the GSF correction to the ISCO frequency [see Eq. ( 85 1] . Parenthetical figures 
show the uncertainty in the last displayed decimals. 



Following Refs. [T()l[55], we parametrize the 0(v) — 0(q) correction to the ISCO frequency by the dimensionless 
parameter Co, such that 

(To! + m 2 )n isco = 6~ 3 / 2 [1 + C a v + 0{v 2 )} , (85) 
where O, lsco is the physical frequency of the ISCO, defined with respect to an "asymptotically flat" time t, and where 
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6 3 / 2 is the dimensionless frequency of the unperturbed (geodesic) ISCO around a Schwarzschild black hole. Reference 
[16] obtained an analytical expression for Cn in terms of the values of a(x) and its first two derivatives at the ISCO: 



Cb = §a(l/6) + l-^, 



where a (1/6) denotes the combination (see next section) 



a(l/6) = a(l/6) + io / (l/6) + io"(l/6). 



(86) 



(87) 



On the other hand, by considering the minimum of the function e(x), Ref. [22j obtained an analytical expression for 
Cn in terms of the values of the first two derivatives of the redshift function zsf(x) at the ISCO: 



Cn 



4^2 



^sf(V6) 



4f(V6) 



It is easily checked that the link (36) (found in Ref. [23]) between a(x) and zsf{x) transforms Eq. (86) into Eq. (88) 



Our goal here is to obtain a more accurate value for Cn based on our new, improved GSF data. Given Eqs. (86) 



and (87), this task amounts to accurately determining a(x) and its first and second derivatives [or, equivalently, the 
perturbation h^ L (x) and its first and second derivatives] at x = g. We comment that in the method referred to 
as second in Refs. [TTJ [TH] (which is modelled upon the scalar-field analysis of Ref. [H]) one essentially also requires 
the second derivative of the metric perturbation from circular orbits, although in this method one derivative is taken 
with respect to the field point (to construct the GSF) and only the second is taken with respect to the orbital 
radius (to consider the effect of a small-eccentricity variation) . From a computational point of view, we hence expect 
both methods (Le Tiec et al.'s [52] and the second method [TTJ [TU]) to perform at a roughly equivalent level. Our 
improved accuracy in Cn will come primarily from using much more accurate numerical data based on the efficient 
frequency-domain code of Ref. [45] . 

As a first attempt, we may simply use our global analytic model(s) for a(x) to read off estimated values for a(l/6), 
a' (1/6) and a"(l/6), and hence for a(l/6) and Cn- For example, our accurate 16-parameter model (#14 in Table [TTT| ) 
and simpler 8-parameter model (54) give, respectively, 



a(l/6) = 0.7958829.. 
a(l/6) = 0.7958860.. 



(using model #14), 
(using model #24), 



(89) 



and 



Cn = 1.2510153. 
Cn = 1.2510199. 



(using model #14), 
(using model #24). 



(90) 



Both values of Cn are consistent with the results of [TTl [19] [22] . However, placing an error bar on our predictions for 
a(l/6) and Cq requires a more careful analysis. Also, it is reasonable to expect that more reliable values for a(l/6) 
and Cn could be extracted from local analysis of the data near x = | rather than relying on global fits. We proceed 
by presenting such a local analysis. 

First, let us give some consideration to the question of the optimal functions for the local analysis near x = g. 
Naively, one may expect either a(x) or cie(x) [or cle{x), or even zsf{x) itself] to be equally suitable for a local fit 
near the ISCO, because all these functions are perfectly regular there. However, one should recall that the rate of 
convergence of the Taylor expansion in x — | about the ISCO (which can be measured by its radius of convergence), 
depends on the global smoothness of the function. For instance, while as (a;) and cle{x) are both bounded and 
continuous (by construction) not only on the closed interval on < x < | but even in larger intervals, say < x < ^, 
the function a(x) itself [as well as a(x) = a(x) / (2x 3 )] blows up at the LR like (1 — Sx)^ 1 ^ 2 . The function zsf(x) blows 
up even faster: like (1 — 3x) _1 . If we assume that the functions we are dealing with can be analytically continued 
in the complex x plane, the radius of convergence of the Taylor expansion around some point xq can (generally) 
be estimated to be equal to the distance separating xq from the nearest singularity, in the complex plane, of the 
considered function [53[ . This suggests that the Taylor expansion (around xq) of all functions having a singularity 
at £ s i ng = | [such as a(x), a(x) or zsf(x)] will have a radius of convergence equal to |a; s i ng — Xq\ = |-| — acol • F° r 
xq = jjr, this yields a radius of convergence equal to |, i.e a Taylor series around the ISCO converging roughly like 
J2n(®( x ~ I))™- By contrast, if we assume, for instance, that the nearest (complex) singularity of the functions a^cc) 
and cie(x) is beyond x s i ng = ^, this suggests that the radius of convergence of the ISCO expansion of these functions 



24 



will be larger than 



= |, corresponding to a series converging roughly like $^„(3(x ~ g))™ (i- e - mucn faster 



than for the LR-singular functions). This suggests that the use of a LR-regular function, such as cle(x), should allow 
for a more accurate determination of the local a(x) derivatives than the direct use of a LR-singular one, such as a(x) 
or z S f{x). 

We have checked this expectation by comparing the relative performances (with respect to local fits) of several 
functions related to a(x). Namely, we considered local fits for the following functions: {di(x):i = 1,...,5} — 
{a, a s , as , a s , cle} , where a s (x) = s(x)a(x) = yl — 3x a(x) and a s {x) = a s (x)f (2a; 3 ). Our analysis began by selecting 
a subset of a,i(x) data (for each i) around x = g. We chose 10 data points on each side of x = |, which, together with 
the x — | point itself, comprise a 21-point subset in the range x € [2/15,0.2], corresponding to r/rri2 £ [5.0,7.5]. We 
least-square fitted each of these di(x) datasets to polynomials in x = (x — g J of degrees varying from 7 to 13, as in 
Eq.(51 1. We judged the quality of each fit by looking at the value Xmin/DoF resulting from each fit. Let Qi(N) denote 
the value of Xmin/DoF corresponding to fitting ai(x) to a polynomial of degree N. In Table|vJ we present the values 
of Qi(N) for i = 1, . . . , 5 (corresponding to the a^s taken in the order of the set defined above) and N = 7, ... ,13. 



N 


ai = a 


a2 = a s 


«3 = ClE 


<i4 = a s 


a 5 = CLE 


7 


1.165 x 10 4 


133.5 


672.9 


2.419 


7.408 


8 


111.1 


2.906 


6.862 


2.188 


2.200 


9 


3.118 


2.387 


2.399 


2.387 


2.388 


10 


2.674 


2.626 


2.629 


2.623 


2.624 


11 


2.831 


2.833 


2.832 


2.847 


2.843 


12 


3.185 


3.186 


3.186 


3.151 


3.160 


13 


3.432 


3.414 


3.432 


3.313 


3.328 



TABLE VI: The values for the "quality of the fit" Qi(N) for the various cn(x) used. N is the degree of the polynomials used 
in the fits. The values in columns 2 through 5 are the Xmm/DoF values for these polynomial fits of degree N. For example, 
<3s(7) is the Xmin/DoF value obtained from a 7 th -degree polynomial fit to cle(x) (with DoF= 21 — 8 = 13). 



The results in Table |V| confirm that the use of LR-regular functions [i.e., a%{x) through a*>{x)] is beneficial in the 
sense that fewer parameters (i.e., lower values of TV) are required to obtain a x 2 /Y)oF of 0(1). Among the LR-regular 
functions, ai(x) = a s (x) stands out as being optimal in that it already reaches x 2 /DoF = 2.419 for N = 7 while for 
this value of TV all the other CLi{x) fare worse, and importantly much worse in the case of the unregularized original 
function a(x), which has x 2 /DoF > 10 4 for N = 7 (and needs at least N = 9 to be considered a good fit). Note also 
that the second-best function is 05(2;) = c\e(x). 

Our strategy, therefore, is to use the function a s (x) for our local analysis at the ISCO. Based on the values presented 
111 Table [Vj above, we use a,4,(x) = a s (x) with N — 8 to compute a(l/6), a'(l/6), a"(l/6), and a"'(l/6) (the latter will 
be needed later) by analytic differentiation of the best-fit model multiplied by 2x 3 j\J\ — 3a; [to translate back from 
a s (x) values to a(x) values]. This yields the following results: 



o(l/6) = 0.0260941094800(93) 
a" (1/6) = 12.00689379(28), 



a' (1/6) = 0.6164354346(12) 
a'" (1/6) = 204.788188(53). 



(91) 



Here the 2-digit error bars refer to the last 2 decimals of each quantity. These errors have been obtained from the 
covariance matrix of the polynomial regression. For instance, the above procedure gives for a'(l/6) an estimate of the 
form Coao + ciad, where &o, cii are the coefficients in the fitting polynomial Pfit(i) = 5o + Sjx 
are coefficients obtained by using the chain rule in differentiating a(x) = 2x 3 a s (x)/\/l 

c§cr o + 2coCict i + cfcn, where 



clnx n , and Co, c\ 



?>x at x = 1/6. This yields a 



squared error on a' (1/6) given by a 2 a , — c§croo + 2coCii7oi + c 2 ci 1 , where an are the elements of the covariance matrix 



coming out of the least-square fit. (Note that here, we are treating the data points as random Gaussian variables 
centered on the values listed in Tables |VIII 



Inserting the values of Eq. (91 1 into Eq. (87) we obtain 



and IX with the variance given by the errors in the tables. 



a (1/6) = 0.795883004(15) [from local fit for cl s (x) 



(92) 



corresponding to [using Eq. (86)] 



Cr 



1.251015464(23) [from local fit for a s (x)}. 



(93) 



Note that the error bar on a(l/6) (and therefore on Cq) is dominated by the error on a"(l/6). As above, the errors 
in these quantities have taken into account the correlations described by the covariance matrix. Actually, we find 
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1/2 — 

that these correlations are rather mild, the largest of which equaling 002/ [000022] — ~0.777 when normalized. As 
a check, we repeated this analysis with our second best LR- regular fitting function, namely qe(x) for N = 8. The 



results agreed with the ones listed in Eq. (91 1 well within the error bars indicated there. For example, the central 
value for a"(l/6) obtained from a local fit to &e(x) is 12.00689372. We also used the non-LR-regular function a(x) 
to repeat the above local analysis using now N = 10 (see Table |v| for why). As expected, it led to larger errors, but 
the central values it gave agreed with the ones listed above within the (larger) error bars entailed by the use of a(x). 
For example, the central value for a"(l/6) obtained from a local fit to a(x) is 12.00689374(54). 

In addition, we also used 7- and 9-point stencil (midpoint) methods, applied to the functions di(x) defined above 
to extract the same derivatives independently. The results for a(l/6) and its derivatives from these stencil methods 
also agreed with the above results within the error estimates corresponding to the stencil methods [which happen 
to be substantially larger, especially as the order of the derivative increases, than those given by the local fits to 
a s (x) or cle(x)]. For example, the 7-point stencil method applied directly to the unrenormalized function a(x) gives 
a"(l/6) = 12.006890(7), where the error was computed considering that the stencil method estimates the derivatives 
as a weighted sum of data points, each of which is, as before, treated as a Gaussian random variable with the variance 
equal to the error listed in Tables |VIII| and |IX| 



As a partially independent check on the value of Cq, we repeated the above analysis working with the variable 



zsf(x). We first constructed a dataset for zsf{x) using Eq. (12 
(with N between 9 and 11) with 7- and 9-point stencils appi 
4^(1/6). We obtained 



with (34)], and then combined a local-fit procedure 



ied to zsf(x) to estimate the values 2^ (1/6) and 



z'(l/6) = 3.10379963(1), z"(l/6) = 22.056551(6), (94) 
where the errors were now estimated from the dispersion between the local fits and the stencil estimates. Hence, using 



Eq. |88|) yields 

Cq = 1.2510155(4) [from local fit and stencils for zsf(x)]. (95) 

Summarizing, all of the above results are consistent with each other within their own errors. We a priori consider 
that it is likely that the most accurate results are those obtained from using the first method above i.e. local fits 
to the LR-regular a s (x). Indeed, by looking at the difference between the data points and the fits for a s (x), one 
sees that they fluctuate in sign in a quasi-random manner across the entire data set. This suggests that the local 
fit is an effective way of averaging out these fluctuations over the 21 data points around the ISCO. However, the 
minimum x 2 /DoF for the local fit is somewhat above unity (namely 2.188 for the best fit used above). [Similarly, 
for our preferred global fit the minimum x 2 /DoF was 4.77]. This hints that, as already mentioned, our error bars on 
the data points have been somewhat underestimated. To be on the conservative side, we simply suggest that all our 
errors bars be uniformly doubled. In particular, this means that we recommend using as our preferred final results for 
a (1/6) and Cq the following values: 

a(l/6) = 0.795883004(30), (96) 



C n = 1.251015464(46). 



(97) 



Our final result (97), which in rounded numbers reads 1.25101546(5), is fully consistent with the value currently 
available in the literature, and it adds 4 significant digits to it. See Table |V| for a comparison. 



VI. ON DETERMINING THE EOB POTENTIAL d(x) 

In this section we discuss the determination of the 0{y) piece d{u) of the second EOB potential D(u]v) 1 defined 
through 

D{u;u) = l + vd{u) + 0{v 2 ). (98) 

[Here, as usual, we use d(u) to denote a functional form.] We present numerical results for d(u) on u < |, i.e., outside 
the ISCO as well as on the ISCO itself, where a certain subtlety occurs. These results are obtained from a combination 
of numerical data and analytic fits. We then comment on the extension of the function d(u) beyond the ISCO. Our 
discussion extends upon and improves the similar discussion presented in Ref. [2"3"] . 



2G 



Reference [TB] obtained a relation involving d(u), the function a(u) (and its first and second derivatives), and the 
function p(u) describing the 0(v) precession effect in slightly eccentric orbits (at the circular-orbit limit). The function 
p{u) is defined for stable circular orbits through 



1- 6a; + vp{x) + 0(u 2 ), 



(99) 



where Q r is the t-frequency of radial oscillations about the circular motion, and 57 is the usual azimuthal i-frequency 
of the circular orbit. As discussed in [TBJ, the definition of p(x) can be extended to include unstable circular orbits 
(i.e., to the entire regime < x < g where timelike circular orbits exist) by replacing the squared radial frequency 
f2j! with (minus) the appropriate squared Lyapunov exponent associated with the growth rate of perturbations of the 
unstable orbit. The said relation between the functions d(u), a(u) and p(u) is given by [16] 



p(u) = Au I 1 



1 - 2u 



a(u) + (1 — 6u)d(u), 



(100) 



where we have introduced [consistent with Eq. (87) above] the shorthand notation 



a(w) = a(u) + ua'(u) + 2 U (1 — 2u)a"(u). 



(101) 



This relation is valid over the entire range where timelike circular orbits exist, i.e. for < u < |. 

The function p(u) was computed numerically in Ref. [41] for u < | , and an analytic fit for it over the corresponding 
domain was obtained in Refs. [141 1181. This, in conjunction with an analytic fit for a(u), allows one to obtai n the 
function d(u) on < u < g via Eq. (100). Reference [53] proposed computing d{u) simply through solving Eq. (100) 
with respect to d(u): 



d(u) 



1 



1 - 6m 



p(u) — a (it) + 4m 



1 - 2m 



1 



(102) 



Note, however, that this expression is formally singular at the ISCO, where 1 — 6m = 0. This singularity should in 
principle be removable (as also discussed in [23J), but the presence of the divergent factor (1 — 6m) _1 makes it difficult 



to evaluate d(u) numerically in the immediate neighbourhood of the ISCO. Thus, the expression (102), as it stands, 



is in practice ill-suited for describing the behavior of d(u) across the ISCO (where this function is expected to be 
perfectly smooth). 

To overcome this difficulty we can use an ISCO-local analysis, as we have don e in the preceding section. An 
expression for d(l/6) can be obtained simply by evaluating the M-derivative of Eq. (100) at u = g [or, equivalently, 
by using de rHopital's rule in Eq. (102)], assuming d(u) is smooth across the ISCO. One finds 



d(l/6) 



p'(l/6)-a'(l/6) 



8V2 



(103) 



which was first derived in Ref. [TB]. Note that a' (1/6) involves all first three derivatives of a(u) at u = i: 



a'(l/6) = 2a'(l/6) + -a"(l/6) + ^'"(1/6) 



(104) 



Hence, calculating d(l/6) requires knowledge of a'(l/6), a"(l/6), a'"(l/6) and p'(l/6). Highly accurate values for the 
former three quantities were given in Eq. (91) above. These values give a'(l/6) = 16.612290(3). 



To obtain the value p'(l/6) we fit to a sample of p(x) data points just outside the ISCO, as was done in Ref. [TS] 
but using a denser sample of data for better accuracy. We obtain 



p'{l/6) = 12.70(5). 



(105) 



This is consistent with the estimated value of 12.66 quoted in Ref. [TH]. [We prefer here to compute p'(l/6) based 
on an ISCO-local analysis, rather than extract this value from any of the global fit models of Refs. [HI [H], since the 



local analysis is likely to produce a more accurate result.] With these values, Eq. (103) gives 



d(l/6) = 0.690(8). 



(106) 
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The uncertainty in this result is entirely dominated by the numerical error in p'(l/6); unfortunately, the latter value 
comes from eccentric-orbit GSF calculations, which are of a relatively limited accuracy. 

Table VII lists all p(x) data available to date, combinin g res ults from Refs. [TB] and [13] • For each data point we 
display t he c orresponding value of d(x), computed via Eq. (102 1 with our accurate 16-parameter model for a{x) (#14 
in Table III I. At the ISCO itself we quote the value obtained above. Figure [8] shows a plot of the d(x) data. Note 
that the error bars on d{x) are predominantly due to the numerical error in the p data, which is larger than the model 
error in our a{x) fit [itself based on very high accuracy a(x) data]. For that reason, it is "safe" to use our analytic fit 
formula for a(x) rather than the actual a(x) data, as we have opted for here. 
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TABLE VII: Available numeric al da ta for p(x) (collected from |14l 1181 l4l] h and quasi-numerical values for the O(v) EOB 
function d(x ), obtained from Eq. ( 102 1 based on the p data in conjunction with our accurate 16-parameter model for a(x) (#14 



in Table 



IIII. 



fCn, using the accurate value we 



The value of p at the ISCO was determined from the simple relation p(l/6) 
have obtained for Cn in Eq. (971 above. The value of d at the ISCO is quoted from Eq. (106 1. Note that the r values in this 
table are exact, while the x values are computed as the inverse of these exact values. 



The data in Table 
methods of Ref.fTHl. 



VII 
We" 



can be used as a basis for an analytic fit model for d(u) over < u < g, e.g., using the 
leave this to future work; we expect that the new, frequency-domain GSF method of Refs. 
[451 149) could soon provide much more accurate data for p(u), that will enable a more reliable and accurate fitting. 
It is of importance to extend the computation of p(u) beyond the ISCO, in order to facilitate the computation of 
d[u > 1/6). This should be possible in principle based on the existing GSF computational framework, although the 
details are yet to be worked out and implemented. 

Such developments would allow one, in particular, to study the behavior of d(u) at the LR. It is interesting to 
speculate about this behav ior b ased on the form of Eq. ( 102 1 and what we already know about the LR-behavior of 
a(u). From Eqs. {47]) and ( |l0l| ) we have a(tt) - §(1 - 3m)" 5 / 2 as u — > | (with £ 

5/2 



§(1 - 3u)~ 5/2 as u -> § (with C ~ 1)- 
we must have that d{u) diverges as d(u 



Hence, if the divergence of 



p(u) at the LR is weaker than oc (1 — 3u) 

will indicate below an argument suggesting that d(u) indeed has a strong divergence oc (1 — 3u) -5 / 2 



1/3) - |(1 -3u)- 5 / 2 . We 
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FIG. 8: The 0(y) EOB function d(x). We plot here the quasi-numerical data from Table VII normalized by the leading-order 
PN term 6x 2 for convenience. 



VII. LIGHT-RING BEHAVIOR AND IMPLICATIONS FOR EOB THEORY 

A striking result of our sub-ISCO GSF computation is the finding that, as one approaches the LR, i.e. as x — >• | 
or u — > |, h^ L (x) blows up proportionally to (1 — 3a;) -3 / 2 , while, correspondingly, a(u) has an inverse square-root 
singularity: a(u) w 0.25(1 — 3m) -1 / 2 . This (apparent) singular behavior a priori raises an issue concerning the domain 
of validity of the EOB formalism, or that of the GSF expansion. In previous EOB work, it was always tacitly assumed 
that all the functions parametrizing the EOB Hamiltonian, i.e. A(u;v), D(u;v) and Q(u,p v ,p r ;is), were smooth 
functions of u = M/VeoB; from the weak-field value u — up to, at least, the Schwarzschild horizon value u = |. A 
smooth behavior of the EOB radial potentials across u = | seems a priori necessary for allowing the EOB formalism 
to describe, for instance, the head-on (or near head-on) coalescence of a (large mass ratio) binary down to Teob ~ 2Af. 
Is the observed singular behavior a(u) « 0.25(1 — 3u) -1 / 2 a signal that something pathological happens in the EOB 
formalism around the radius teob ~ 3M, or is it an artefact of formally trusting the first [O^ 1 )] term in the GSF 
expansion beyond its physical domain of validity? Let us start discussing this issue by considering the physical origin 
of the LR singularities in h^ L (x) and a(u). 



A. Physical origin of the light-ring divergences of h^ L (x) and a(u) 

We first explain the singular behavior h^' L (x) oc (1 — 3x)~ 3 / 2 in terms of heuristic technical considerations. [The 
weaker singular behaviour a[u) oc (1 — 3it) -1 ' 2 then follows from the structure of Eq. (42).] Note that we are interested 



in the regularized field (x ), which is obtained from the full, retarded perturbation h^ u (x ) by subtracting the 
Detweiler- Whiting S-field h^(x x ) (and evaluating the result on the particle's circular orbit). We will consider in 
turn the LR behavior of the full and singular fields, and argue that the singular behavior of h^ L (x) is inherited from 
the full field (while the S-field remains bounded at the LR). 

Consider first the full Lorenz-gauge metric perturbation /i^„(x A ), which is sourced by the stress-energy tensor 
of particle 1, i.e. T^ u (x x ) = mi(-g)~ 1/2 J dru iJ -u v 5{x x - y x {r)) = {~g)~ 1/2 J miu^dy" 5(x x - y x {r)), where r, 
y x {r) and u x = dy x /dr are, respectively, the particle's proper time, trajectory and four- velocity. As the particle 
approaches the LR (along a sequence of circular orbits), the (non- vanishing) components of its 4- velocity (in any 
frame at rest with respect to the background Schwarzschild frame of m^) tend towards infinity proportionally to 
u° = dt/dr = (1 — 3x) -1 / 2 . Therefore, the (non- vanishing) components of T^ v too will tend to infinity proportionally 
to (1 — 3a;) -1 / 2 . In other words, we can write T^ v = (1 — 3x)~ 1 / 2 T tJ - l/ where all the components of the "renormalized" 
stress-tensor stay bounded as particle 1 tends to the LR. Correspondingly, we can write = (1 — Sa;) -1 / 2 /^, 
where the "renormalized" metric perturbation is sourced by T^ v , so that can be written as the convolution 
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of a suitable tensorial Green function with T M „. The latter convolution might introduce an additional, milder singular 
behavior in the LR limit |54j , but it is unlikely to alter the leading-order power-law blow-up oc (1-3X)- 1 / 2 . Then, the 
value of the redshift-related scalar h^ u , which contains two extra factors (1 — 3a;) -1 / 2 coming from the two contractions 
with the four-velocity (we assume here the "constant" off-worldline extension of the four-velocity discussed in Sec. 
II B I, is expected to blow up near the LR proportionally to (1 — 3a;)~ 3 / 2 



Let us next consider the S-field h^(x x ). Near the particle, the trace- reversed counterpart of this field, h^{x x ) = 
h^(x x ) — \h s - L {x x )g tlu {x x ), has the leading-order form h^(x x ) ps Amiu^u u /e, where e is the invariant orthogonal 
geodesic distance between x x and the worldline. When contracting this local expression with u^u v , the two factors 
of u M disappear (by virtue of u^u^ = — 1) and we are left with h^u^u" ~ 4mi/e. Now, e scales proportionally to 
m m near the LR, so that h^, and hence also h^u , finally scales proportionally to the inverse of u M near the LR, i.e., 
huu (1 — 3x) +1 / 2 . Note that the above "cancellation" of factors (due to u^u^ = — 1) does not occur in the case 
of the full field h^ u (x x ), which is obtained by a global integral that generally does not yield a result proportional to 
u»u v . 

In summary, we find that the LR behavior of the regularized difference field h^ L (x x ) = h^ u (x x ) — hf^(x x ) is 
dominated by the LR behavior of h^ u (x x ), i.e., it is naturally expected to blow up near the LR proportio nally 



to (1 — 3x)~ 3 / 2 . We can say that the blow ups h^ L oc (1 — 3a;)~ 3 / 2 , and correlatively [according to Eq. (42 1] 
a(u) oc (1 — 3u) -1 / 2 , as particle 1 tends to the LR are simply rooted in the corresponding power-law blow-up of the 
components of the 4-velocity near the LR: oc (1 — 3a;) -1 / 2 . 

Having understood this simple technical origin of the LR behavior, we can reformulate it in a physically more 
transparent way. Instead of parametrizing (as is usually done) the strength of first-order GSF effects by means of 
the rest-mass mi, one can say that the source of the perturbation is better measured by the conserved energy of 
the small mass, say E\. Indeed, we recall that £\ is given by a hypersurface integral of the contraction of the stress- 
energy tensor with the time-translation Killing vector k^d/dx^ — d/dx° i.e. E\ = — J k^T^dS^. We note also that 
£i = — mig^k^u" , clearly exhibiting the fact that Ei measures the eventual growth of the components of u M . In the 
case of circular orbits, it takes the simple form E\ = miE c i TC (uo) where uo = mijrQ is the background gravitational 
potential at the considered orbital radius r , and where we used the notation 

, x 1 - 2u , 

Scire (u) = -==. (107) 
VI — oil 



For later conceptual clarity we have here added a subscript "circ" to the function E(u), already defined in Eq. (48 1. 

We conclude that a physically transparent way of interpreting the blow-up of a (it) near the LR is to say that the first- 
order GSF correction a(u) oc m,i/m 2 actually grows proportionally to the ratio of the conserved energy E\ — miE c i IC (u) 
of the small mass to the large mass m-2 : i.e. a(u) ~ qE c - uc (u). We have already used this reformulation above to 
factor out the singular function E c i lc (u) from the GSF-data-derived function a(u). As for h^ L , its stronger blow-up 
near the LR is equivalent to saying that h^ L (x) ~ z~ 1 a(x) ~ q z^ 1 E c i IC (x) , with z = 1 — 3x as usual. 



B. Physical domain of validity of GSF results. 

Before discussing the impact on the EOB formalism of the LR-divergent behaviors of h^ L (x) and a(u), let us 



address the issue of the physical domain of validity of the GSF approximation. We already mentioned (in Sec. II C ) 
the evident condition that a first-order [O^ 1 )] GSF calculation makes sense only if the conserved energy of the small 
mass, £i, is parametrically smaller than that of the large mass, £2 ~ nri2- hr the context of circular orbits, this leads 
to the necessary condition 

Ei = mi E cilc f 7 ^) < m 2 , (108) 



where E c i IC (u) is the function defined in Eq. (1071. 



Actually, this necessary condition is not sufficient for the consistency of the GSF expansion. Indeed, though we do 
not yet have a second-order GSF calculation of h^ L and a(u), one can physically estimate that second-order GSF 
effects will (at least approximately) modify the zeroth-order (geodesic) expression Ei — mi-E C irc(^) used in the 
condition ( 108[ ) above for the energy of mi by including the back reaction of mi on the background metric. Therefore, 



one expects that a more accurate version of the above necessary condition will roughly read 

mi Scire ( m2 ^, Q C£l ) < m 2> ( 109 ) 
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where we have modified the zeroth-order gravitational potential uq — m 2 /ro by replacing m 2 by m 2 + c£i, where c is 
some constant of order unity. However, if we now look at the crucial square-root contained in the singular denominator 



of condition (109), it reads 



3^_ 3 3 = «/l-3^-c^A. (HO) 
r r V r r m 2 



Near the LR, the hopefully more accurate condition ( 109 1 (which approximately takes into account second-order GSF 
effects) is consistent with the first-order GSF condition ( 108| ) only if we have 



i.e. 



£ L« Zo ^l_ 3 m2 , (111) 
m 2 r 



v « (zo) 3/2 - (H2) 



This is much stronger (near the LR, i.e. when z — > 0) than the condition (108) which corresponded to qE c [ IC (uo) <C 1, 
or q w v -C J~Zq. We can also note that the approximate inclusion (following the pattern used above) of third-order, 
and higher-order, GSF effects do not a priori seem to require stronger constraints on v. Indeed, our treatment above 
essentially consisted of considering a first-order fractional modification of the "effective background mass" (say near 
the LR, where r7i 2 /r ~ 1 does not introduce an independent small parameter) of the type m 2 — > m 2 [l + Ci^-E c ; rc (uo)]. 
We can generalize this treatment by considering higher-order GSF contributions (proportional to higher powers of 
Ejjmz) leading to a replacement of the type m 2 —¥ m 2 [l + ci^i? c i rc (uo)+c 2 (^i?circ(ito)) 2 + - ' ' ]■ However, the condition 



( 111 I is such that all the higher-order terms c n (vE^ TC (uo)) n , (n > 2) are consistently smaller than the first-order one 

Cl^circ(uo)- 

Let us also remark in passing that another way to understand the necessity of the stronger consistency condition 



(111) is to notice that it is tantamount to requiring 



h^ L ^Zo 1 —«l. (113) 

This makes sense because h^ t L yields the first-order GSF perturbation of the proper time of particle 1. More precisely, 
the regularized proper time of particle 1 reads 



Clearly, it makes sense to expand the squareroot in powers of q only if the actual magnitude of h^ L is small compared 
to 1. 

Summarizing so far: first-order GSF results near the LR are a priori physically meaningful only in a limit where 
the ratio i>/(zo) 3 ^ 2 tends to zero. 

C. Singular light-ring behavior as a coordinate singularity in the EOB phase space 

The findings of the previous subsection imply that the physical implications of the mathematically divergent LR 
behavior a(u) - (1 - 3u)~ 1/2 of the 0{v) piece of the EOB radial potential A(u, v) = 1 — 2u + ua(u) + 0(v 2 ) 
are less dramatic than they seem to be at first sight. Indeed, in the domain where we can trust the derivation 



of this result, i.e. under the condition (112), the first-order GSF contribution to A(u, v) remains small, namely 
va(u) ~ vjyfz « z « 1. In add ition , the first-order GSF contribution to the it-derivative of A(u,v) also remains, 
under the consistency condition (112), much smaller than unity, namely ua'iu) ~ v/z 3 / 2 <C 1. On the other hand, 
if we consider the second u-derivative of A(u,u), it will be of the form A"(u,v) = va" (u) + 0{v 2 ) - v/z b / 2 , which 



increases so much near the LR that the condition ( 112 ) is compatible with arbitrarily large values of A"(u, v) as one 
approaches the LR. 

Therefore, it is a priori possible that the near-LR behavior a(u) ~ (1 — 3u) -1 / 2 only corresponds to some mild 
type of physical singularity at the LR, and that higher-order effects in v will smooth out the shape of the A potential 
into some sort of "boundary layer" near the LR. However, the appearance of a formal square-root singularity makes 
it unclear how the function a(u) can be extended across the LR to radial arguments u > 1/3. Indeed, the formal 
analytic continuation of (1 — 3U)" 1 / 2 leads to imaginary values of a(u) when a > 1/3. 
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However, independently of the possible role of higher-order effects, we are faced with the mathematical fact that 
the 0{v) piece in A(u, v), i.e. the value of the ^-derivative dA(u, v)jdv at v = has a singularity a(u) ~ (1 — 3u)~ 1//2 . 
Is this inescapable mathematical singularity signalling the presence of some real singularity in the EOB formalism? 
We think instead that it does not correspond to any physical singularity of the EOB dynamics, but is simply a 
coordinate singularity in phase-space, which can be avoided by a suitable (symplectic) phase-space transformation. 
Indeed, a somewhat similar, formally singular LR behavior of the EOB A(u) potential has already appeared in another 
EOB work, in a problem where one can see how this apparent singularity can be avoided by a suitable phase-space 
transformation. We are alluding here to a recent analysis by Bini, Damour and Faye |50j of tidal effects in comparable- 
mass binary systems, based on an effective action approach, completed by an EOB reformulation. Reference |50| found 
that the perturbative description of these effects leads, within the standard EOB description of circular orbits, to a 
radial potential of the form 

A(u, v, mt) = A 2pp (u, v) + [i T a T (u, v) + O(^), (114) 

where the term /it o-t(u,u) denotes the additional contribution, coming from tidal interactions, to the "two-point- 
particle" EOB radial potential Ai m {u,v') . Here [It symbolically denotes a generic tidal parameter (actually there 
is a sum over a set of tidal parameters: [It = /J.^, /if \ 1^3 > ' ' ' )• ^ was proven in Sec. VI of Ref. |50j that, in the 
extreme-mass-ratio limit v« 1, the (various) tidal contributions a?{u, v) are singular at the LR. More precisely, they 
formally blow-up as ot{u,v — 0) ~ (1 — 3u)~ 1 when u — > |. Let us sketch how this singularity in ax(u) can be 
avoided by a suitable phase-space transformation that replaces it with an alternative regular description. 



First, it should be noted that when comparing Eq. ( 114 ) and the GSF-expanded result A(u. v) = 1 — 2u + v a(u) + 
0{v 2 ) the analogy is between a perturbation expansion in powers of \xt ("tidal expansion") and a perturbation 
expansion in powers of v (GSF expansion). To make the argument more crisp, let us actually set v to zero in Eq. 
( 114[ ), i.e. consider tidal effects on a body of mass ra\ <C ni2- [It is in this limit that the results of Ref. [50] which we 



shall use below could be rigorously proven.] Let us now recall how Ref. [50 derived the presence of a LR singularity 
in a^(u). This was done in essentially two steps: (i) computation of the additional effective action due to tidal effects, 
and of its Hamiltonian formulation; (ii) reformulation of this original Hamiltonian perturbation as a contribution to 
the standard EOB Hamiltonian in the limit of circular motions. The result of the first step is that tidal effects add 
to the squared effective Hamiltonian (H c g(u,p v ,p r )) 2 a new contribution which, for general orbits, is quartic in the 
(effective) momenta p^, say 

5 T H 2 s (u,p v ,p r ) = pTq{u,Pp) + 0(a4), ( 115 ) 

with 

q(u, Pfl ) = C^ x (u) PtiPu p K p x , (116) 

where the tensor C is a smooth function of u (in particular it is regular at u = [In the above expressions the time 
component po is meant to be replaced by (minus) the unperturbed effective Hamiltonian.] For instance, for the leading- 
order quadrupolar tidal effects the tensor C pvK \ is proportional to a symmetrized version of (1 — 2u) R a q v R KaX , 
where R a jlj ^ u is the background curvature tensor. 

It is the second step taken in Ref. [50] (namely the reformulation into the standard EOB Hamiltonian form) which 
actually introduced the singular behavior at u = |. To explain this point, let us first recall that the standard form 
of the EOB Hamiltonian, introduced in Ref. [J], consists of imposing some specific restrictions on the momentum 
dependence of the squared effective EOB Hamitonian. Namely, it should have the form 

Hc S {u,p v ,p r ) = A(u;v) (^l+p^u 2 + B ^.^ +Qr CS tr{u,p ip ,p r ;v)j , (H7) 

with the specific form p^u 2 of the term quadratic in p v , and with the restriction that the mass-shell deformation 
term Qrestr vanish quartically in the limit of small radial momentum: 

Qr CS u(u,p v ,p r ->• 0) = 0(p4) . (118) 

Reference 0] showed, at the 3PN accuracy, how such a restricted form can be reached by applying a suitable sym- 
plectic transformation of EOB phase-space variables (q*,Pi) — > (q' z ,p'i)- At the 3PN accuracy, it was found that 
Qiestr(u,p,p,Pr) did not depend upon p v , and was given by Q3PN(u,p r ) — 2(4 — 3v)vu 2 pf- We, however, expect that 
Q restr will involve a dependence on p v at higher PN orders (see the discussions in [16 and the appendix of [23] ). 
Let us also recall that this standard EOB form is well tuned to the description of near-circular orbits. For instance, 
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as discussed in [IB], it allows one to describe the dynamics of small-eccentricity orbits only in terms of the radial 
potentials A(u) and B(u). 

Reference |50j used the fact that, for the special case of circular motions, the value of the original first-order Hamil- 
tonian perturbation ( |115[ ) must coincide with the corresponding value in the transformed phase-space coordinates, 
(q',p ! ). In accord with the results we shall derive below in a GSF context, this led to the following link between the 
original, LR-regular momentum-dependent perturbation (116) and the momentum- independent tidal perturbation of 
the standard radial A potential, entering (114): 



MtOU„ 



l 



pi u 2 



(119) 



Here, we have added a prime on the radial variable appearing on the left-hand side as a reminder that the function 
ar(u) belongs to the transformed phase-space coordinates (q' ,p'). However, at the first order in fix we are considering, 
v! can (and actually shoul d) b e identified with the dummy variable u used on the right-hand side (RHS). Of crucial 
importance in the result ( |119 ) is the fact that all the quantities on the RHS should be evaluated along the one- 
parameter sequence of circular motions. This means that the momentum components p^ on the RHS are to be 
replaced as follows: ~po is replaced by E c [ TC (u), Eq. (107), while the tangential component of the spatial momentum, 
say p lt = p v u is replaced by 



pf %u) = upf c {u) 



1 - 3u 



(120) 



Here, and in the following, we shall often use, without explicating the change of notation, scaled EOB variables, such 
as Pi = pf hys //i or r = r phys /M. We also remove the label EOB from the radial coordi nate r for notational simplicity. 



Finally, as the (transformed) inverse radial variable u tends to | , we see that the link ( 119 ) implies that the potential 
ot(u) blows up proportionally to 



0(PJ) 



(pf c (u)) 2 ~(l-3 U y 



(121) 



The main point we wanted to emphasize by explicating the appearance of the pole singularity in ar(u) (as u — > |) 
found in Ref. |50j . was that its origin was the growth, for large values of the tangential momenta p n , of the original, 
momentum-dependent Hamiltonian perturbation ( 115[ ) which was, to start with, perfectly regular near (and across) 
U = |. As we shall see in more technical detail below (in the GSF case), it is actually the change of phase-space 
variables (q,p) — > (q' ,p'), needed to go to the restricted, standard EOB Hamiltonian form, Eq. ( 117), that is responsible 
for introducing a singularity in ot(u). 

Let us now apply the same reasoning to our GSF context, i.e. using v as a perturbation parameter, instead of \it 
in the tidal case above. To clarify this application (as well as what was at work in the tidal example above) we shall 
also show how to explicitly construct the phase-space transformation, say T : (q l ,pi) — > (q'EOBiP'f B )j needed to go 
to the restricted, standard EOB Hamiltonian form ( |117[ ), parametrized by two radial potentials A(u) and B(u), and 
an EOB mass-shell deformation function Q r estr(u,p v ,p r ) constrained to be 0(pf) when p r — > 0. In the discussion of 
Rcf. -TOl recalled above, the transformation T was implicitly used, but its explicit form did not enter. 

Let us start from the unperturbed (squared effective) EOB Hamiltonian 



H C «o( U ,Pv'Pr) = A o(u) 



1 



B (u) 



(122) 



with Aq(u) = 1 — 2u, E>o(u) = (1 — 2u) _1 , and consider that, in some original phase-space coordinates (q l ,Pi) (say, 
coordinates directly related to a Lorenz-gauge calculation) GSF effects modify it by adding a new contribution of the 
general form 



S u H 2 s (u,p v ,p r ) = vq(u,p v ,p r ) + 0(y 2 ). 



(123) 



We assume that, in the original phase-space coordinates (q l ,Pi) = (r,ip,p r ,p v ) (for planar motions), the perturbed 
Hamiltonian in Eq. (123) is a smooth function of {q l ,Pi), and, in particular, that no divergence occurs when u crosses 



the value |. On the other hand, we allow for a general (unrestricted) dependence of q{u,p v ,p r ) on the momenta. In 



other words, the sum H^ so (u,p v ,p r ) + S^H^ s (u,p v ,p r ) is not assumed to be of the standard EOB form, Eq. (117). 

Let us now look for a symplectic phase-space transformation, say T : (q l ,Pi) — > {q'EOBjP'f° B )> which simplifies the 
form of the original squared effective Hamiltonian, H^ SQ (u,p v ,p r ) + S^H^ s (u,p v ,p r ), by putting it in the standard 
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EOB form of the type displayed in Eq. (117), with the restriction ( 118 ) . We shall sketch here how one can formally 
construct the needed symplectic transformation T within our GSF-perturbative context. (In previous work, T was 
constructed within a PN-perturbative context.) T can be taken to arise from a (perturbed) generating junction of 
the form g(r,p v ,p r ) — h , p r e r (r,p 2 ,p 2 ) + 0{v 2 ). [In preparation for some technical developments below, we use the 
space-time symmetries of the two-body dynamics to infer that the time- reversal-invariant, scalar quantity e r can be 
written as a function of p 2 = (q x p) 2 and p 2 = (q • p) 2 /r 2 .] The 0(v) phase-space transformation generated by 
g(r,p v ,p r ) is effected via a Poisson bracket, namely 8 g X = {X, g}, on any phase-space function X(q,p). From a 
technical point of view, let us formally consider that all the dynamical functions of interest are expanded in powers of 
p r , while keeping the dependence on r and p v exact. (This differs from the usual PN-perturbative construction of T 
which essentially uses a multiple expansion in powers of p r ,Pu = up v and u.) A simple calculation (using {r, p r } = 1, 
{PnPip} = 0, etc.) then shows that 5 g p r — {p r ,g} — 0{p r ). In other words, the transformation T respects each 
order in the expansion in powers of p r , namely <5 9 p™ = 0(p"). For simplicity, we shall focus here on the terms of 
zcroth-order in p r . These terms are already quite nontrivial. It can be verified that the reasoning indicated below can 
be straightforwardly extended to higher powers of p r . 

At zeroth-order in p r , the original Hamiltonian perturbation v q(u,p v ,p r = 0) + 0(v 2 ) depends on two independent 
phase-space variables, namely u and p^. Our aim here is to show how a suitable generating function g can reduce 



the general dependence of q(u,p v ,p r = 0) on u and p v to the special one entering Eq. (117). At zeroth order in 
p r , the only relevant changes in phase-space variables are that in r and p v . [The change in ip is irrelevant as the 
relevant dynamical observables do not explicitly depend on (p.] A simplifying feature is that 8 g p v — {p^,g} = 
(because dg/dcp = 0). Let us then consider the change in r: 5 g r = {r,g} = dg/dp r . This is easily found to be 
S g r = v e r (r,p 2 j ,p 2 ) + 0(p 2 ) + 0{v 2 ). In other words, to zeroth order in p r , and to first-order in v [i.e., modulo 
corrections of 0{p 2 ) + 0{v 2 )] we have 5 g r — ve r {r,p 2 v , 0). Note that the change in radial coordinate is more general 
than a simple (configuration-space) coordinate transformation Sr — C{ r ) m that 5 g r depends on both r and p v . This 
phase-space dependence of S g r is crucial for allowing the transformation T to reduce the original contribution Eq. 



(123) to the standard canonical EOB form. For convenience, we shall work in the following with the corresponding 
change in u = 1/r, i.e. 5 g u — —S g r/r 2 , and denote it as S Q u — ve u (r,p 2 n ) [mo dulo corrections of 0(p 2 ) + Q(v 2 x] 



The condition on g is that it transforms the sum of Eq. (122) and Eq. (123) into the standard form Eq. (117), with 



some modified potentials A(u) — Aq( u) + va{u) + 0(v 2 ), B{u) = Bq(u) + vb(u) + 0(v 2 ) : and some restricted 0{v) 



mass-shell term Q res tr satisfying Eq. (118). Written explicitly, this condition means that S g H^ [ u,p v , p r ) must be 



equal to the difference between (123) and a GSF perturbation of the standard EOB Hamiltonian (117). The latter 



GSF perturbation has the structure 

^cff standard^, ?V) = M») U + P^ + J^) "&(«) + ^ 

where the contribution 0(p^) comes from the restricted 0{v) mass-shell term Q rcs t r - At lowest order in p r , and after 
division by the condition on g reads (with j = p v to ease the notation) 

U g [(1 - 2«)(1 + fu 2 )} = q(u,j,p r = 0) - o(«) (1+M (125) 



If we introduce the short-hand notations 

, 2 2u(l - 3u) 

e = 11,-2 2 e ( r >r 

and 

u{u,f) = — - — , (126) 



: 2s - ff( M >J>r = ) 

l + j 2 u 2 



the latter condition explicitly reads (after dividing both sides by 1 + j 2 u 2 ) 

e(u,j 2 )(j 2 ~3 2 ilc (u)) = a(uj 2 ) - a(u), (127) 

where 

j ^ U) S u^u) (128) 

denotes the function of u which describes the value of j = p v along the sequence of circular orbits. [One must carefully 
distinguish the general, independent phase-space variable j from the specific function j c i rc (u).] 
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From the condition (127) on g, one first deduces [by taking the limit j — > j 2 lTC (u) on both sides], that the (a 



priori unkown) value of the canonical perturbed EOB potential a(u) corresponding to the original non-canonical 



perturbation in Eq. ( 123 ) is given by 

a{u) = a(u,j cilc (u)) = -r—^2 T\~t - 12 9) 



Then, we derive the value of the function e r (r,p 2 ,p 2 = 0) determining the transformation g (at lowest order in p r ) 
from 

2e"(r,j 2 ) = a _ a(u, j 2 ) - a{u, j 2 irc (u)) 

Note that the last expression on the RHS defines, despite the appearance of a denominator that vanishes along the 
phase-space curve j 2 — j 2 irc (it), a smooth function of u and j 2 along this seemingly singular curve. This follows from 
the fact that a is a smooth function of its second argument. [Here, we are using the fact that, if f(x) is a smooth 
function of x, g{x,y) = {fix) — f(y))/(x — y) is a smooth function of the two variables x and y, even in the vicinity 
of the diagonal x — y.] 



The result in Eq. ( 129 I is the GSF-perturbation analog of the (tidal-perturbation) result of Ref. [50] cited in Eq. 



(1191 above. It is this result which explains the appearance of LR singularities in the "standard" EOB potential 



a(u) when starting from a LR-regular o rigin al non-standard perturbed EOB Hamiltonian, Eq. (123). Indeed, if the 
original 0(v) GSF perturbation in Eq. (123) is regular in phase-space (including near u = |), but grows as pf when 



the components pi of the momenta get large, we see from Eq. (129) that such a growth at large momenta in the 
original phase space will lead, after the transformation T, to a corresponding growth of the purely radial function 
a (it) as u — > | on its (transformed) u axis of the type 

a («) « i lh c{ 2 ] ( \ - (jcUu)r- 2 ~ , L 2)/2 - ( i31 ) 

1 + W JcircW (l — 3u)( n 2 >' Z 

Reciprocally, if we reason backwards, our construction above of the 0{u) generating function g can be used (as 
we shall explicitly discuss in the next subsection), when starting from the singular standard 0{v) EOB potential 
a(u) ~ (1 — 3m) -1 / 2 , to transform it away, and to replace it by a regular, (unrestricted) momentum-dependent 0{v) 
contribution to the EOB Q function. 

We therefore conclude that our finding above of a LR singularity in the perturbed standard 0(v) EOB potential a(u) 
prob ably 55j originates from an everywhere-regular unrestricted perturbed 0{v) effective Hamiltonian v q(u,p v ,p r ), 



Eq. (123), which grows cubically (i.e. cx p^) as p v — > oo. Note that the new, transformed u axis corresponds to the 
original phase-space variable v! = u + ve u (r,p 2 ) + 0(v 2 ). Note also that the result Eq. (130) formally determining 



the transformation T to the canonical EOB form becomes ill-defined as u tends to | . Worse, if we approximate [for 
large j 2 and large j 2 irc (u)} q(u,j 2 ) as ~ j 3 , so that a(u,j 2 ) = q(u,j, 0)/(l + j 2 u 2 ) ~ j, we see from Eq. (130) that 



■2 -2 

u -1-1 J ~ Jcirc J Jcirc /ioo\ 

e ~ J Jcirc - 2 _ ,2 ~ • , • ■ 
J Jcirc J < Jcirc 

Not only is this result blowing up as either j c i rc or j gets large, it is actually only well-defined above the LR, i.e. when 
r > 3 or u < |, because it contains the square root j c i IC (u) = (u(l — 3w)) -1 / 2 . Therefore, the transformation T, and 
the corresponding standard EOB potential a(u), are (probably) only defined when u < | . 

One can also check that the reasoning above can be extended to higher orders in p r . In particular, at order 0{p 2 ) it 
determines the value of the second "standard" EOB potential 8 u B(u) = B{u;v)—Bq{u) = vb(u)+0{v 2 ). A preliminary 
study of the 0{u) contribution to the standard B potential indicates that it blows up like b(u) cx (1 — 3u)~ 5 / 2 as 
U —> As a consequence, d(u) will have the same type of divergence near the LR. 

We can summarize our conclusion by an analogy. For many years, researchers in general relativity have been 
mystified by what they called the "Schwarzschild singularity", namely the fact that the standard Schwarzschild 
metric is singular at r = 2M, notably because g rr = (1 — 2A//r)~ 1 blows up there and then changes sign. It 
was only in the 1960s, notably through the work of Kruskal, that it became clearly understood that this "r = 2M 
Schwarzschild singularity" is a singularity of the standard Schwarzschild coordinates, which can be gauged away by a 
suitable transformation of the spacetime coordinates, including a necessary mixing of space and time coordinates. Our 
conclusion is that the singularity a(u) ~ (1 — 3u) -1 / 2 we found is, somewhat similarly, only due to the a singularity 
of the standard phase-space coordinates used in the EOB formalism. This "phase-space-coordinate singularity" can 
be gauged-away by a suitable symplectic transformation, necessarily mixing coordinates and momenta. 
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D. Impact of our findings on the EOB formalism 



What conclusions should we draw from our GSF sub-ISCO results for the EOB formalism? One possible reaction 
would be to modify the standard EOB strategy that concentrates all the near-circular dynamical information into the 
effective metric [parametrized by the two radial potentials A(u) = A (u) +va(u) + 0(v 2 ) and B(u) = B (u) + ub(u) + 
0(V 2 )], and to allow the third EOB function, Q(u,p r ,p v ; v), to participate in the description of near-circular orbits, 
by relaxing the constraint that Q(u,p r ,p v ; v) vanishes when p r — > 0. (In the analogy with the "r = 2M Schwarzschild 
singularity" case, this reaction is the analogue of modifying the standard Schwarzschild coordinates by allowing the 
use of a more general spacetime gauge fixing.) Let us sketch how this could be done. For this purpose, it is convenient 
to introduce a special notation for the piece of the squared effective Hamiltonian contributed by the Q function. Let 
us denote 



Q(u,p Vl p r ;is) = A(u;v)Q(u,p Vl p r ;v) 
so that we have a simple linear decomposition of the squared effective Hamiltonian: 



H cS {u,P v ,Pr) = A{u;v) (l+p v u ) + 



A(u; v) 
B(u; v) 



pl + Q(u,p v ,p r ;is). 



(133) 



(134) 



Note that we are no longer adding a subscript "restr" to Q or Q. Indeed, we are no longer imposing the constraint 
(118 1, but we allow a general momentum dependence in Q and Q. 



Let us now consider the GSF expansion of the squared effective Hamiltonian in Eq. (134), corresponding to the GSF 
expansion s of A , B and Q, namely A(u) = A (u) + va{u) + 0(v 2 ) and B(u) = B (u) + vb(u) + 0(v 2 ) and, consistently 
with Eq. ([123]), Q(u,p v ,p r ; v) = vq(u,p v ,p r ) + {v 2 ). It yields Hl a {u,p v ,p r ) = H 2 S0 (u,p v ,p r ) + 8 y Hl^{u,p v ,p r ), 
where the zeroth-order term is given in Eq. ( |122 ) above, and where the 0{v) perturbation is given by 



8 v Hl s (u,p v ,p r ) = va(u) (l+p 2 v u 2 ) + , 



a(u) 



b(u 



A Q (u) 



p 2 r +vq(u,p <p ,p r ) + 0(v 2 ). 



^B (u) " v 'B 2 (u), 

As above, let us focus on the terms in the Hamiltonian which survive in the limit where p r — > 0: 
S^H^iu^p^^r = 0) = i>a(u) (l +p 2 p u 2 ) +vq(u 1 p Vl p r = 0) +0(v 2 ). 



(135) 



(136) 



This contrasts with the corresponding p r —> limit of the perturbation of the standard, restricted EOB Hamiltonian 
which would only contain the first contribution, linked to the perturbation of the standard A potential. 

In previous sections, we (following, in particular, Ref. |23j ) have interpreted the numerical GSF data by assuming 
that we were working within the context of a standard EOB Hamiltonian. It was within this context that we found 
a standard a(u) potential of the form aE{u)E c { rc {u) . In other words, we interpreted the GSF data in terms of the 
following perturbation of the standard-gauge EOB Hamiltonian 



^H 2 Ssi . andald (u,p v ,p r = 0) = v a E (u)E ciTC {u) (l + p%u 2 ) + 0(v 2 ). 



(137) 



where E c i IC (u) is the function of u alone defined in Eq. (107) 



Let us now discuss the many ways in which the latter, LR-singular standard EOB Hamiltonian ( 137 1 can be traded 
off for an everywhere-regular Hamiltonian of the general, non-standard form (136). From Eq. (127), the criterion for 



two Hamiltonians to be equivalent (at zeroth order in p r ) modulo a symplectic transformation is simply that their 
numerical values agree along the sequence of circular motions, i.e. when p v = p™ c (u). This criterion leaves many 
possibilities for transforming (137) into an equivalent, but LR-regular Hamiltonian. 

The simplest way of doing so is to replace the problematic factor E c i IC (u) = H c g (u, p" TC (u) , p r = 0) in (137) by the 
v — > limit of the full (non circularly reduced) effective EOB Hamiltonian, i.e. the square-root of Eq. ( 12"2FTnamely 



H e so(u,p vl p r ) = a A (u) 1+p^u 2 + 



B (u) 



(138) 



This leads, when considering for simplicity the p r — hypersurface in phase space [56] . to the following first non- 
standard possibility for replacing (137): 



K H es( U 'P<P>Pr = °) = V a E(u)H cS0 (u,p v ,p r = 0) (l +p 2 p U 2 ) + 0(h> 2 ). 



(139) 
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Considered as functions over phase space, the two pert urbed Hamiltonians in Eqs. (137) and (139) are very different 



functions. In particular, the new Hamiltonian Eq. ( 139 1 is regular across u = = [57], because it is constructed from the 



"regularized" function Oe(u), extrapolated beyond u = |, as discussed in Sec. III. In addition, 8' v H^ s {u,ptp,p r = 0) 
vanishes at u — \ at least as fast as \J\ — 2u. [Actually, as we saw above that the regularized function qe(u) is likely 



to change sign near u = ^, the new Hamiltonian Eq. (139), being proportional to ob(m)VT 
(1 



2u) 



3/2 



near u 



2it, nearly vanishes as 



The p v dependence of Eq. (139) is oc (1 



Pi u 2 



,3/2 



This is consistent with our general conclusion above that any 
LR-regular form of the 0(y) perturbed Hamiltonian must grow as p 3 for large momenta. However, it was argued in 
Ref. [4] that it is most natural for the Q contribution that it be 0(p 4 ) for small momenta. It is easy to accommodate 
such a requirement by considering new reformulations of the naive possibility Eq. ( 139 ) within the "equivalence class" 
of perturbed Hamiltonians that numerically agree along the sequence of circular motions [i.e. when p r — 0, and 

V% =icirc( U )]-_ 

Let us consider the following phase-space function 



2 

k(u,p v ) = (l-2u) ^ a 

1+PtpU 2 



(140) 



It is easily seen that, along circular orbits [when p 2 v = j 2 i IC (u)] the phase-space function k[u,p v ) is numerically equal 
to 1. O n the other hand, as a function of momenta it is 0{p 2 ) for small momenta. We can then defin e a reformulation 
of Eq. (137) which is both regular at u = | and 0(p 4 ) for small momenta, by multiplying Eq. (139) by the square of 



k(u,p v ), say 



6"H 2 s (u,p v ,p r = 0) = va E (u) (k(u,p v )) H e s (u : p v ,p r = 0) (l + p 2 J u 2 ) + 0{y 2 ) 



giving, explicitly, 



S^H 2 s (u,p v ,p r = 0) = va E {u) u 2 (1 - 2uf' 2 



+ 0(v 2 ). 



(141) 



(142) 



Note that the global p dependence of this new Hamiltonian is quite different from the previous one [namely p A / \Jl + p 2 
instead of (1+p 2 ) 3 ^ 2 , where p stands for p u — up v ], though (by consistency) they both grow like p 3 for large momenta. 
Note also that this new Hamiltonian has a faster vanishing near u = \, namely oc oe(u)(1 — 2u) 5 ^ 2 



More generally, one could use a phase-space transformation which trades off the GSF result Eq. ( 137 ) for a perturbed 
Hamiltonian such that 



f n H* s (u, Ptp ,p r = 0) = v f(u,k(u,p v )) H eS0 (u,p v ,Pr = 0)(l+plu 2 ) +0(u 2 ), 



(143) 



where f(u, k) is any function such that f(u, 1) = oe(u). This clearly leaves a lot of freedom in the definition of such a 
LR-regular version of the standard result (137). Note that the phase-space function k(u,p v ) tends to the finite limit 
(1 — 2u)/u as p v — > oo, independentl y of b eing restricted to the sequence of circular orbits, so that the whole class of 
regular perturbed Hamiltonians Eq. ( 143 ) grows as p 3 for large momenta. 



One should not be surprised by the existence of such a large freedom in the formulation of a regular, non-standard 
Hamiltonian. Indeed, given any specific, LR-regular Hamiltonian, its transform by an arbitrary regular symplectic 
transformation 7~ reg will generate a new, regular Hamiltonian. The fact that, at zeroth order in p r , such a regular 
Treg can introduce an arbitrary function of two variables f(u, k) [constrained only along a certain curve in the (u, k) 
plane] is linked to the presence of the arbitrary function e r (r,p 2 ,p 2 = 0) in the p r — > limit of a generating function 
g(r,Pip,p r ) = v p r e r (r^^jp 2 ) + 0(v 2 ). We leave to future work the detailed generalization of our considerations to 
higher orders in powers of p r , and simply note that, anyway, such a generalization involves [at 0(v)] an arbitrary 
function of three arguments [corresponding to e r (r : p 2 ,p 2 ))- We just recall here that part of the success of the EOB 
formalism consists of finding good ways of trimming down this large gauge freedom to parametrize the dynamics in 
terms of the minimum number of relevant functions. 

We have just explained how to make a full use of our sub-ISCO results, without being restricted by their singular 
behavior at the LR, by relaxing the "standard gauge fixing" of the EOB formalism. However, it should be noted that 
such a modification of the current EOB formalism is really only needed if one wishes to describe the dynamics of 
ultra-relativistic quasi-circular orbits (p — > oo) near u = ^. By contrast, the original motivation for, and main use of, 
the current EOB formalism is to describe the dynamics of mildly-relativistic radiation-reaction-driven quasi-circular 
orbits. Such orbits stay close to the sequence of (stable) circular orbits down to the ISCO (i.e. for < u < h), and 
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then strongly deviate from the sequence of unstable orbits that formally continue to exist when g < u < |. Indeed, 
though the "plunging motion" that follows the radiation-reaction-driven quasi-circular inspiral remains approximately 
circular (i.e. withp 2 <§; p 2 , see Fig. 1 in Ref. [3]), its path in phase-space (q,p) drastically deviates from the phase-space 
location of unstable circular orbits. In particular, the angular momentum p v of a plunging orbit stays approximately 
equal to its value jlf£° = Vl2 + 0(v) when it crossed the ISCO, while the formal adiabatic sequence of p v value s 
along the unstable circular orbits is given by the function j C irc(w) defined as the square root of the RHS of Eq. ( 128| . 
In particular, as one gets near u = | the two phase-space points ((/piungeiPpiunge) an d (qdrcPdic) become infinitely 
far apart. This infinite phase-space separation (in the p v direction) is both the cause (as we have seen above) of 
the divergence of the standard a(u) as u — > |, and an indication that the latter divergence is not of direct physical 
relevance for describing (as the EOB formalims aims to do) the dynamics of plunging orbits. Actually, a proof of the 
capability of the standard 3PN-accurate EOB formalism (as defined in [3]) to accurately describe the dynamics of 
(comparable-mass) coalescing black hole binaries down to the light-ring has been recently given in Ref. |15l . Figure 1 
of the latter reference shows in particular that the (uncalibrated) standard 3PN-EOB prediction for the £(j) curve 
agrees remarkably well with the NR one down to the radial location of the LR (indicated as the leftmost vertical line 
in the figure). By contrast, the £ (j) curve corresponding to the formal adiabatic sequence of circular orbits starts to 
exhibit a strong, and increasing deviation from £ NR (j) after the crossing of the ISCO (see the dash-dotted line in the 
latter figure). 

In view of this effectiveness of the standardly gauge-fixed EOB formalism for the description of the dynamics of 
mildly-relativistic binary systems, it might be useful to set up a minimal way of modifying the EOB formalism so as 
to incorporate our GSF sub-ISCO results. [In our above analogy, this is like continuing to use standard (or nearly 
standard) Schwarzschild coordinates when describing a system for which the coordinate singularity at r — 2M is not 
interfering with the physics one is interested in.] Above we have discussed ways of entirely trading off the 0(v) piece of 
the radial A potential for an equivalent momentum-dependent Q-type contribution. Actually, the latter momentum- 
dependent Q-type contribution (growing cx p^ for large momenta) is only needed for describing ultra-relativistic 
motions, while the standard A-type contribution is a simple and effective way of describing mildly-relativistic motions 
above the ISCO. One can then conceive of a mixed scheme, where the dynamics is described partly by a certain 
radial v a$(u) potential, and partly by a vq contribution, with the vao{u) potential playing the leading role during 
the inspiral, and the v q contribution taking over only during the (late) plunge. For instance, considering as above 



the terms remaining when p r — > 0, if we constrain the vq contribution to have the same p-dependence as in (142 1 



namely p 4 / \/l + p 2 , we can use 

CE*s(u,p v ,Pr = 0) = va (u) (1+pIu 2 ) + 

4 

''' (a E (u) u 2 (1 - 2uf' 2 - a Q (u) u 2 (1 - 2u) 3 / 2 (l - 3u) 1/2 ) + 0{y 2 ). (144) 



Here one can choose ao(it), which corresponds to a {/-deformed radial potential A(u; v) = Aq(u) + v oo(m), at will. 



For instance, one could choose ao(u) = 2 it , so that the p / yl + p 2 -type contribution in Eq. (144 1 starts, when 



u is small, proportionally to u , i.e. at the 3PN level [by contrast to Eq. (142 1 which starts like u p^j + u 2 p 2 , 

which corresponds to the 2PN level]. Alternatively, one could choose an ao(u) which stays very close to the "exact" 
standard one (a^(w)i? C i rc (u)) up to some value u — uq, and then deviates from it when u > uq, and stays regular 
across the LR. Such a choice would ensure that the p A / 1 \fl + p 2 -type contribution in Eq. (144) has a negligible effect 



when u < uq, and starts modifying the dynamics only when u > uq. For instance, one could choose a value of Uq 
between the ISCO and the LR. This would allow one to make full use of our new strong-field results on a(u) up to 
u = Uq < 3, essentially without modifying the EOB formalism up to u = uq. [It is with this program in mind — of 
defining some simple, accurate ao(u) approximation to a standard (w) — that we have given (above) accurate estimates 
for the first three derivatives of a standard (u) at the ISCO.] The LR-regularized p 4 j ' \J\ +p 2 -type contribution in Eq. 



(144 1 would then only affect the end of the plunge which follows the crossing of the ISCO. 

We leave to future work a study of the performances and relative merits of the various possible completions of the 
EOB formalism discussed above, as well as a discussion of the needed extra terms of order 0(p 2 ) (S-type contributions) 
and 0(p 4 ) (old, standard Qrcstr-type contributions). In this respect, we note that it would be very valuable to be able 
to use GSF data on plunging motions to directly extract EOB-useful information about the plunge dynamics taking 
place after the crossing of the ISCO. Alas, the current state of development of GSF theory (namely the lack of explicit 
0(y 2 ) results) does not allow one to extract gauge-invariant information from the calculation of the gauge- variant self- 
force along a plunging orbit. In a related vein, we note that, even if we go back to the case of exactly circular orbits, 
our current GSF calculation of the first-order only, 0(v), contribution to the (standard) radial A potential is quite 
insufficient for allowing one to construct an estimate of the function A(u; v) able to accurately describe the dynamics 
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of comparable-mass binary systems. Our current best-bet knowledge of the function A(u\ v) for comparable-mass 
(non-spinning) systems (i.e. v ~ 4) has been obtained by: (i) introducing [9] a two-parameter family of putative 
A-functions incorporating current analytical (and GSF) knowledge, and then (ii) best-fitting, for each available value 
of u, the corresponding EOB-predicted waveform to NR waveforms [21 HD] . The results of these EOB/NR fits indicate 
that the function A(u; v) cannot be accurately described with only a linear dependence on v, i.e. a function of the 
form A(u; v) = 1 — 2u + va(u). The latter fact was explicitly discussed in the Conclusions of Ref. [16], especially 
around Eq. (8.2) there. Even if one disregards the indications from such EOB/NR fits about the need for 0{v 2 ) terms 
in A(u, v), our GSF results on the behavior of a(u) near u — | give another hint about the need for and importance of 
such terms. Let us just quote two illustrative examples which involve mildly strong-field effects |58j in the dynamics 
of a small mass around a large mass. 

First, there is the issue of the existence of an ISCO. We know that it exists in the geodesic limit v — > 0, and is 
located at r phys = 6M, i.e. u = g. We should a priori expect that such a mildly strong-field phenomenon continues 
to exist (when neglecting radiation-reaction effects) as the symmetric mass ratio increases to values of order 4 . In the 
(standardly gauge-fixed) EOB formalism, the condition for the existence of an ISCO (defined as the condition for the 
existence of an inflection point in the effective potential describing the radial motion) is |16] 

A(ui S co,^) = 0, (145) 

where misco is the looked- for inverse radius, and where the function A (it, i>) is defined as 

A(u,v) ee 2A(u,v)A'{u,v) + 4u(A'(ii,i/)) 2 - 2u A(u, v) A"(u, v), (146) 

with the prime denoting d/du. By mathematical continuity with the solution uisco — g which exists when v — 0, 
the condition (1451 will certainly admit a solution of the form Uisco = \ + 0(v) in a neighborhood of v — 0. We 



can now explore the physical need for terms of order 0(y ) in the function A(u; v) by studying the range of values 



of v where a solution of ( 145 1 continues to exist under the assumption that the A potential is assumed to be exactly 
linear in v, i.e. given by the formula A hn (u; v) = 1 — 2u + v ogsf(u) with our above-determined sub-ISCO a function, 
o-gsf{u) — aE{u) E c i rc (u). We find that such a solution exists only for a rather small neighborhood (0, i/ max ) of 
v = with v max 0.108. Note in passing that a radial A potential of the type A lm (u;^) = 1 — 2u + voqsf(u) 
also has the unpleasant physical consequence of predicting the existence of some stable rest position, at some radius 
r, with fixed values of 9 and (p. Indeed, the effective radial potential H 2 S (u,p v ; v), considered as a function of 
the radial coordinate u — 1/r, for any fixed angular momentum p v , reduces in the case where p v = to simply 
A lm (u]v) = 1 — 2u + ^aGgp(u), which exhibits, for any non-zero value of is, a minimum at some value of u < | 

(because <2gsf( u ) +°o as u —> | ). 

Let us give a second illustrative example for the undesired physical consequences of keeping only the term linear in 
v in the A potential. It concerns another mildly strong-field effect. In the geodesic limit, one of the unstable circular 
orbits plays a somewhat preferred role. It is the one located at u = | (i.e. r phys = AM). It has a (specific) angular 
momentum p v — 4, and a vanishing binding energy, i.e. E = 1. This marginally bound orbit is the end point of the 
special zero-binding zoom-whirl orbit which starts, in the infinite past, with zero kinetic energy at infinity (but with 
the non-zero angular momentum p v — 4) and ends up, in the infinite future, "whirling indefinitely" around the large 
mass. As in the case of the ISCO, one would a priori expect that such a mildly strong-field phenomenon will continue 
to exist (when neglecting radiation-reaction effects) as the symmetric mass ratio increases to values of order | . The 
condition for such a zero-binding circular orbit to exist has been written down (in the EOB formalism) in Ref. |16j . 
It reads 

Z(u*,v)=0, (147) 
where it* is the looked-for radius, and where the function Z(u, v) is defined as 

Z{u,v) ee A{u,v) + ~uA'(u,u) - A 2 (u,v) . (148) 



We know that the condition (147) admits the solution u* — \ when v = 0. By mathematical continuity, it will 



certainly admit a solution of the form = | + 0(y) in a neighborhood of v = 0. Like in the case of the ISCO, we 



can now explore the range of values of v where a solution of ( 147) continues to exist under the assumption that the A 



potential is assumed to be exactly linear in v. We find that this leads to a much stronger constraint on the magnitude 



of v than in the case of the ISCO. Namely, we find that a solution of (147) exists only for a very small neighborhood 
(0, v' max ) of v = with f^ ax ~ 0.035. For larger values of v the growth of oqsf(u) in the interval | < u < | prevents 
the continued existence of a zero-binding circular orbit (as well as of the corresponding zero-binding zoom-whirl orbit). 
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These striking physical consequences of neglecting 0(v 2 ) terms in A(u, v) suggest that the 0{v 2 ) contribution to 
A(u, v) is also divergent near u = |, but has a negative sign. In addition, one expects it to diverge proportionally to 
(1 — 3«)~ 2 so that the ratio between the 0(v 2 ) contribution and the 0{v x ) one is of order — 3u) 3 / 2 , as suggested 



by the previously derived consistency condition (112 1. We leave to future work a detailed study of the higher-order 
GSF contributions to A(u, v). 

VIII. CONCLUSIONS AND OUTLOOK 

We have computed the conservative piece of the gravitational self-force (GSF) acting on a particle of mass m x as 
it moves along any (stable or unstable) circular geodesic orbit around a Schwarzschild black hole of mass mi ra\. 
Our main results and conclusions are as follows. 

(1) We numerically computed the function h^ L (x) = h^i L u^u v , where h^ L (oc mi) is the regularized metric 
perturbation in the Lorenz gauge, v) 1 is the four- velocity of m\ in the background Schwarzschild metric of m2, and 
x = [ Gc~ 3 (mi + m2)fl] 2 / 3 is a dimensionless measure of the orbital frequency Q. Our results are collected in Tables 
VIIl| and VIII in Appendix |a") The fractional accuracy of our numerical results ranges between 10~ 10 and 10~ 8 for 



most data points, and never gets worse than 10 -5 (except at a single point, closest to the LR). Our results improve 
on previous calculations both in accuracy (except for very small values of x, x < 1/200), and in range. In particular, 
our work is the first to explore the unstable orbits between x ~ \ [slightly below the innermost stable circular orbit 
(ISCO) located at x = |], and the light ring (LR), located at x = g. 

(2) We particularly studied the behavior of h^ L (x) just outside the LR at 1 = | (i.e., r = 3Gm 2 /c 2 ), where the 
circular orbit becomes null. We found that h^ L (x) blows up like h^ L (x |) ~ _ |^C(1 _ 3x) -3 / 2 , where 

We argued that the divergence of h^ L (x — > ^) can be understood from the divergent behavior of (some of) the 
components of the four-velocity u M near the LR. 

(3) Using a recently discovered link [23] between h^ L (x) and the piece a(u), linear in the symmetric mass ratio 
v = m\mij(jn\ +ra 2 ) 2 , of the main radial potential A(u, v) — l — 2u + is a(u) + 0(is 2 ) of the effective one body (EOB) 
formalism, we computed from our GSF data the EOB function a(u) over the entire domain < u < ^. Our results for 
the function a(u) improve on previous calculations both in accuracy, and in range. In particular, our work is the first 
to explore the behavior of the EOB potential a(u) as u — > |. We found that a(u) diverges like a(u) ~ jC(l — 3u) -1 / 2 

(where £ ss 1) at the light-ring limit, u — > Q) . 

(4) We then considered the energy-rescaled function oe(u) = a(u)/E(u), where E(u) = (1 — 2u)/\/l — 3u is the 
(specific) relativistic energy of mj in the background Schwarzschild black hole of mass m 2 3> mi. This energy-rescaled 
function has a finite limit as u — > ^, but seems to have a weak singularity ~ cq + (1 — 3u)(c 1 1 ° s In 1 1 — 3u\ + ci) there. 
We gave several high-accuracy global analytical representations of ob(m) that incorporate all the presently known 
post-Newtonian (PN) analytical information about it, and essentially reproduce all our numerical results within their 
numerical errors. We think that our analytical models of oe(u) give a reasonably accurate representation of the 
behavior of that function even beyond the range (0 < u < |) where GSF data can compute it, say in the range 
I < u ^ \- See notably the curves for models 13, 14 and 19 in Fig. [5J which show the Newton-rescaled function 
o-e{u) = a#(u)/2w 3 . In other words, we think that our GSF calculations give us, for the first time, valuable information 
about the truly strong-field regime u = GM/c 2 r < |. 

(5) Using our accurate analytical fits of the EOB potential a(u), we computed global analytical representations of 
the 0(v) pieces in the functions giving the total energy and total angular momentum of a binary system in terms of 
the frequency parameter x. We found that these 0(v) functions have rather strong (negative) divergences near the 
LR, namely ~ — cv{l — 3a:) -2 (with positive constants c). 

(6) The GSF-induced, 0{v), shift in the value of the orbital frequency of the innermost stable circular orbit (ISCO) 
has been a touchstone for comparing various analytical descriptions of binary dynamics. Using a multi-pronged 
analysis of our accurate new data, we have been able to improve the computation of the 0{v) ISCO shift by four 
orders of magnitude — see Table [V] We have also expressed our improved result in terms of the combination a (1/6), Eq. 



(87), of derivatives of the EOB potential, thereby providing a direct, accurate way of calibrating the EOB formalism in 
the v — >• limit — see Eq. ( [96] ). In addition, for further helping the construction of 0(^)-accurate EOB Hamiltonians, 
we have giv en accurate numerical estimates of the values of a(u) and its first three derivatives at the ISCO point 



U=\: Eq. (|9T1) 



(7) In previous work we used GSF data on slightly eccentric orbits to compute a certain linear combination of a(u) 

and its first two derivatives, involving also the 0[y) piece of a second EOB radial potential D(u) = 1 + v d(u) + 0(v 2 ). 

Combining these results with our new accurate global analytic representation of a{u), we numerically computed d(u) 

on the interval < u < \. 

— 6 
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(8) Our finding of an inverse-square-root singularity a(u) oc (1 — Su)^ 1 ^ 2 in the 0(v) EOB potential seems to 
put in question the domain of validity of the GSF expansion, and/or that of the EOB formalism. We addressed 
both issues in detail. First, we argued that 0{v) GSF results are physically reliable near the LR only in a limit 
where the ratio ~ 3it) 3 / 2 tends to zero. This limit allows for the unboundedness, near u = |, of second (and 
higher) u-derivatives of the EOB potential A(u, u), and thereby signals the presence of some type of singularity of the 
(standard) EOB formalism at u = | . However, we argued that the (mathematical) singularity a(u) oc (1 — 3u) -1 / 2 
we found is only a spurious singularity, due to the use, in the current, standard EOB formalism, of some specific way 
of fixing the phase-space gauge freedom. [The a(u) oc (1 — 3u)~ 1//2 singularity is a phase-space analog of, e.g., the 
9rr = (1 — 2M/r) _1 "Schwarzschild coordinate singularity" at r = 2M.] We explicitly showed (at lowest order in p r ) 
how to "gauge-away" the singularity a(u) oc (1 — 3u) -1 / 2 by relaxing the standard, phase-space gauge-fixing conditions 
of the EOB Hamiltonian (namely by allowing the third, Q EOB potential to grow oc p^ when p r = and p v — > oo). In 
addition, we exhibited minimal ways of modifying the current EOB gauge-fixing, which are appropriate when dealing 
with radiation-driven inspiralling and coalescing binaries, rather than with highly unbound ultra-relativistic circular 
orbits near the LR. In order to globally construct these modifications of the current EOB formalism, it is essential 
to make use of our finding that, after factoring E(u) out of a(u), one ends up with a function that is continuous at 
U = ^, and can be naturally extended to larger values of u. 

Finally, let us make the following remarks about some of the future research directions that suggest themselves to 
complete our results. 

(a) First, it would be useful to improve the accuracy of the GSF data near the LR, in order to allow a better 
characterization of the behavior there. This would likely require a reformulation of the mode-sum scheme to achieve 
a more rapid convergence of the multipole mode-sum near the LR, where the standard high-Z behavior is no longer 
applicable. 

(b) It would be interesting to extend the GSF computation of the precession of slightly eccentric orbits so as to 
extend the range of determination of the second EOB potential, d(u), from the current range < u < | to the full 
range < u < | where it is, in principle, computable. 

(c) Several aspects of our work have emphasized the need for an understanding of higher-order terms in the GSF 
expansion, notably terms 0(v 2 ). This provides a motivation for pushing more effort in this direction. 

(d) Our work has also emphasized the importance of being able to extract gauge-invariant dynamical information 
about plunging orbits. 

(e) Finally, it would be interesting to study the performances and relative merits of the various non-standard EOB 
schemes whose necessity is suggested by our work. 
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Appendix A: Numerical data 



We tabulate here the complete set of GSF numerical data used in our analysis. We show numerical values for 
the Lorenz-gauge quantities h^ L (x) sampled at (generally) equal intervals in x in the range 1/150 < x < 1/3, 
corresponding to 3m 2 < r < 150m2- For practical reasons we split the data between Table VIII (x < 1/6) and Table 
IX (x > 1/6). The tables also show the corresponding values of the 0{v) EOB potential a(x) derived from h^ L (x) 
using Eq. (42 1. Our data for x > 1/5 is new. A subset of the x < 1/5 data already appeared in the literature [2T1 [35] 
but at significantly lower accuracy. Reference [24], which was concerned primarily with the PN domain, presented 
a sample of very high accuracy data for the weak-field range 1/500 < x < 1/200. As we are interested here in the 
global behavior (in particular in the strong-field domain), we do not include these high-accuracy large-r points in our 
sample to avoid statistical bias in our \ 2 analysis. All of our data points are consistent with previously published 
results within the respective error bars (where quoted). [To see the agreement with Refs. [21] or [24] one needs to use 
our Eq. (151 in order to convert between the Lorenz-gauge values given in our tables and the flat-gauge values given 
in those sources.] 

To allow for a meaningful x 2 analysis in the present work, it was important for us to obtain a reliable estimate of 
the numerical error in the data points. Our methods for error estimation are described in detail in Refs. [19] and [45) . 
For most data points the error is by far dominated by the uncertainty in the value of the analytically fitted large-Z 
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tail contribution to the regularized metric. Essentially, we use some of the our large-/ numerical data points to fit 
a power-law model using several plausible models (varying over the number of power-law terms and the number of 
data points used for the fit), and use the variance of the results as a rough measure of the tail-fit error. See [H?l 143] 
for more details. We expect this procedure to give us the actual error to within a factor ~ 2 or so. (This is indeed 
confirmed by our \ 2 analysis: we find that the value of x 2 /DoF settles at around 3-4 and does not reduce any further 
upon adding model parameters.) 

The parenthetical figures in Tables [VIII| and |IX| correspond to the error estimates coming from the above procedure. 
For instance, 0.02693868484(1) stands for - 0.02693868484 ± lCT 11 . In this example, ±10~ n describes our "best 
guess" for the numerical error, although a more conservative approach (taking into account the uncertainty in the 
error itself) would perhaps set this at ±2 x 10~ n . Note that the tables also quote (in the fifth column) more "precise" 
values for the numerical errors, given to 3 places right of the decimal point. Strictly speaking, this extra information 
is, of course, meaningless, given the factor ~ 2 uncertainty in the errors; we present it here only for the purpose of 
allowing interested readers to fully reproduce our \ 2 analysis. 
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TABLE VIII: Numerical data (part I). Each row displays data for a circular geodesic with a particular Schwarzschild radius r, 
given in the third column. The first and second columns show the corresponding values of x = m^jr and z = 1 — 3m2/V. [The 
relation x — 1112/r holds only at O(is ), but here we may ignore higher-order corrections because x and r are used as arguments 
(independent variables) for quantities which are already O(v), namely h^ L (x) and va(x).\ The fourth and last columns give, 
resp ectiv ely, our numerical results for the Lorenz-gauge quantity [see Eq. jfij)] and for the 0{v) EOB potential a(x) [see 

Eq. ( 42 1] . In the fifth column we give estimates of the absolute numerical errors in the h^ 4 data. These error values are 
the ones used in our \ 2 analysis, and we give them here in full (showing several insignificant digits) in order to allow readers 
to reproduce this analysis accurately. Our actual error estimates for the individual data points (which we expect to be only 
accurate to within a factor 2 or so) are expressed in the form of parenthetical figures in the fourth and last columns, showing 
our best estimate of the uncertainty in the last displayed decimals. 
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300a; 



r/m 2 



Num. err. 



a(x) 



50.0 
51.0 
52.0 
53.0 
54.0 
55.0 
56.0 
57.0 
58.0 
59.0 
60.0 
61.0 
62.0 
63.0 
64.0 
65.0 
66.0 
67.0 
68.0 
69.0 
70.0 
72.0 
74.0 
76.0 
78.0 
80.0 
82.0 
84.0 
86.0 
88.0 
90.0 
91.0 
92.0 
93.0 
94.0 
95.0 
96.0 
97.0 
97.5 
98.0 
1575/16 
98.5 
99.0 



0.50 
0.49 
0.48 
0.47 
0.46 
0.45 
0.44 
0.43 
0.42 
0.41 
0.40 
0.39 
0.38 
0.37 
0.36 
0.35 
0.34 
0.33 
0.32 
0.31 
0.30 
0.28 
0.26 
0.24 
0.22 
0.20 
0.18 
0.16 
0.14 
0.12 
0.10 
0.09 
0.08 
0.07 
0.06 
0.05 
0.04 
0.03 
0.025 
0.02 
1/64 
0.015 
0.01 



6.000000 
5.882353 
5.769231 
5.660377 
5.555556 
5.454545 
5.357143 
5.263158 
5.172414 
5.084746 
5.000000 
4.918033 
4.838710 
4.761905 
4.687500 
4.615385 
4.545455 
4.477612 
4.411765 
4.347826 
4.285714 
4.166667 
4.054054 
3.947368 
3.846154 
3.750000 
3.658537 
3.571429 
3.488372 
3.409091 
3.333333 
3.296703 
3.260870 
3.225806 
3.191489 
3.157895 
3.125000 
3.092784 
3.076923 
3.061224 
3.047619 
3.045685 
3.030303 



-1.0471854796(1) 
-1.08659952251(6) 
-1.12774434980(7) 
-1.17075624628(7) 
-1.21578528730(6) 
-1.26299706191(6) 
-1.31257466418(8) 
-1.36472098296(8) 
-1.4196613648(1) 
-1.47764670375(8) 
-1.5389570493(1) 
-1.60390583535(9) 
-1.6728448488(1) 
-1.7461701015(1) 
-1.8243287924(1) 
-1.9078276091(1) 
-1.9972426612(1) 
-2.0932314430(1) 
-2.1965473024(1) 
-2.3080570549(1) 
-2.4287625435(1) 
-2.7026089937(1) 
-3.0299861727(1) 
-3.4275063529(1) 
-3.9189791853(2) 
-4.5395895523(2) 
-5.3432356852(6) 
-6.416105450(1) 
-7.903439422(2) 
-10.066672801(2) 
-13.4187934749(9) 
-15.849341504(3) 
-19.093469318(8) 
-23.580398048(7) 
-30.07746738(2) 
-40.08190371(5) 
-56.8876661(5) 
-89.13862(2) 
-118.32926(2) 
-167.13828(2) 
-244.5136(1) 
-260.3517(2) 
-484.6(5) 



1.054 x 10~ 1U 
6.245 x 10" 11 
7.185 x 10' 11 
7.489 x 10' 11 
5.720 x 10" 11 
6.359 x 10' 11 
7.999 x 10' 11 
7.670 x 10" 11 
9.639 x 10' 11 
7.896 x 10" 11 

1.137 x 10" 10 
9.043 x 10' 11 
1.087 x 10" 10 
1.107 x 10" 10 
1.355 x 10~ 10 
1.351 x 10" 10 
1.363 x 10 -10 
1.371 x 10~ 10 
1.373 x 10" 10 
1.375 x 10~ 10 
1.401 x 10~ 10 
1.445 x 10" 10 
1.496 x 10" 10 
1.471 x 10~ 10 
1.627 x 10" 10 
1.609 x 10" 10 
6.309 x 10~ 10 
1.404 x 10~ 9 
1.905 x 10~ 9 
1.559 x 10~ 9 
9.284 x 10" 10 
2.644 x 10~ 9 
8.388 x 10~ 9 
7.359 x 10~ 9 
2.200 x 10~ 8 
4.853 x 10~ 8 
4.832 x 10~ 7 
1.684 x 10~ 5 
1.775 x 10~ 5 

2.138 x 10~ 5 
1.271 x 10~ 4 
1.914 x 10~ 4 
5.398 x 10" 1 



0.02609410950(3) 
0.02821688301(2) 
0.03048093197(2) 
0.03289458866(2) 
0.03546673669(1) 
0.03820686140(1) 
0.04112510844(2) 
0.04423234741(2) 
0.04754024627(2) 
0.05106135426(2) 
0.05480919704(2) 
0.05879838596(2) 
0.06304474249(2) 
0.06756544251(2) 
0.07237918263(2) 
0.07750637433(2) 
0.08296936903(2) 
0.08879272321(2) 
0.09500350907(2) 
0.10163168283(2) 
0.10871052135(2) 
0.12437313326(2) 
0.14234657312(2) 
0.16308580174(2) 
0.18718609087(2) 
0.21544503763(2) 
0.24896018744(6) 
0.2892884360(1) 
0.3387190693(1) 
0.40077307330(9) 
0.48120301414(5) 
0.5312203677(1) 
0.5902619091(3) 
0.6612773504(3) 
0.7488226643(7) 
0.860429954(1) 
1.00975332(1) 
1.2250733(3) 
1.3763418(2) 
1.5789875(2) 
1.828231(1) 
1.872213(1) 
2.357(3) 



TABLE IX: Numerical data (part II), covering the sub-ISCO range 1/6 < x < 1/3. The table is structured in the same way as 
Table Em] 



